© The Authors 2025

1. Introduction

Modern galaxy redshift surveys, such as the Sloan Digital Sky Survey (SDSS; York et al. 2000), Two Degree Field Galaxy Redshift Survey (2dFGRS; Colless et al. 2001; Cole et al. 2005), WiggleZ Dark Energy Survey (Drinkwater et al. 2010), Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013, 2016; Alam et al. 2021), Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016), and Euclid (Laureijs et al. 2011), have enabled us to study the large-scale structure of the Universe, which is mainly traced by galaxies. To fully exploit the next generation of observational data, we need accurate models of galaxy clustering.

In typical galaxy surveys, measurements of clustering statistics, such as the two-point correlation function (2PCF), are more precise at megaparsec scales, which is below the valid scales of state-of-the-art perturbation theory approaches. The fully nonlinear treatment, N-body simulations, is essential for accurately resolving the nonlinear dynamics in structure formation. However, simulations are computationally demanding, posing a significant bottleneck in the exploration of cosmological parameter space.

In contemporary galaxy formation and evolution theories, galaxies form and evolve in dark matter halos. After the distribution of halos is obtained from simulations, we adopt the framework of halo occupation distribution (HOD; Peacock & Smith 2000; Seljak 2000; Berlind & Weinberg 2002; Zheng et al. 2005; Guo et al. 2014) as the galaxy-halo connection model. HOD models have been widely applied to interpret galaxy clustering data (e.g., Zheng et al. 2009; Guo et al. 2015; Zentner et al. 2019). In its basic form, the HOD quantifies the likelihood, P(N|M), that a dark matter halo of mass M contains N galaxies. Together with models of the spatial and velocity distributions of galaxies inside halos, the HOD framework can interpret the clustering of galaxies down to nonlinear scales.

Given the computational requirements of N-body simulations and advancements in machine learning techniques in recent years, emulators have emerged as an increasingly popular tool for modeling galaxy clustering statistics within the HOD framework (e.g., Kwan et al. 2015; Kobayashi et al. 2022; Nishimichi et al. 2019; Zhai et al. 2019, 2023; Kwan et al. 2023; Cuesta-Lazaro et al. 2023, 2024; Burger et al. 2024; Storey-Fisher et al. 2024). The core idea of emulators is to evenly sample the parameter space and apply novel interpolation schemes to create high-accuracy estimates of clustering statistics at any point within that space. We have built a model based on the ABACUSSUMMIT simulations (Maksimova et al. 2021).

The radial positions of galaxies are usually derived from the observed galaxy redshifts assuming a cosmological model. However, these derived distances differ from the true positions due to peculiar velocities, an effect known as redshift-space distortions (RSDs). A traditional way to probe the real-space galaxy clustering and mitigate the RSD effect is through the projected 2PCF, w_p(r_⊥), obtained by integrating the three-dimensional 2PCF, ξ^S(r_⊥, r_∥), along the line-of-sight (LOS; r_∥) direction. In this work, we limit our attention to the real-space galaxy 2PCF, ξ^R(r), building emulators for it and demonstrating that the emulators can recover unbiased cosmological constraints in a likelihood analysis. The modeling of redshift-space clustering will be presented in subsequent studies, where the emulators for ξ^R(r) will serve as one of the model ingredients.

The structure of this paper is as follows: In Sect. 2 we introduce the theoretical background for the galaxy clustering. In Sect. 3 we give a brief overview of the suite of cosmological N-body simulations used in this work and outline how we constructed galaxy catalogs from them. In Sect. 4 we present our neural network emulator. Then, we present our method for recovering cosmological parameters in Sect. 5, before presenting the results in Sect. 6. Finally, we discuss our results in Sect. 7 before concluding in Sect. 8.

2. Theoretical background

The spatial distribution of galaxies can be quantified in terms of the 2PCF, ξ_gg(r), which measures the excess probability of finding a pair of galaxies at a given separation, compared to that of a random distribution. In cosmological surveys, the comoving distance to a galaxy, r(z), is inferred from its redshift via

$$\begin{array}{c}\hfill r(z)={\displaystyle {\int }_0^z}\frac{\mathrm{d}z\prime }{H(z\prime )},\end{array}$$(1)

where H(z) is the Hubble factor. We used natural units where the speed of light is c = 1.

In addition to the Doppler shift resulting from the Universe’s expansion, light emitted from galaxies experiences an additional shift due to peculiar velocities. These peculiar velocities induce anisotropies in the measured spatial distribution of galaxies, known as RSDs. If a galaxy pair is separated by a vector r in real space, it will be separated by a vector s in redshift space. We adopted the distant-observer approximation, where galaxies are assumed to be sufficiently distant such that the angle they subtend with respect to the LOS direction, denoted by $\widehat{\mathit{\ell }}$, is negligible. Under this assumption, galaxies’ peculiar velocity only induces a redshift in the $\widehat{\mathit{\ell }}$ direction. To the linear order in velocity, |v|≪1, the redshift-space separation of galaxies due to peculiar velocities is given as

$$\begin{array}{c}\hfill \mathit{s}=\mathit{r}+\frac{v_{\Vert }}{a(z)H(z)}\widehat{\mathit{\ell }},\end{array}$$(2)

where $v_{\Vert }\equiv \mathit{v}·\widehat{\mathit{\ell }}$ is the peculiar velocity of the galaxy along the LOS, a(z) is the scale factor and H(z) the Hubble factor.

We decomposed the separation vectors r and s into a component parallel and perpendicular to the LOS, s_||,  r_|| and s_⊥,  r_⊥, respectively. Under the distant-observer approximation, RSD effects are isolated into s_||, meaning we have r_⊥ = s_⊥. Consequently, projecting galaxies onto a plane orthogonal to r_|| results in the same distribution of galaxies in real and redshift space, known as the projected correlation function,

$$\begin{array}{c}\hfill w_p(r_{\perp })\equiv {\displaystyle {\int }_{-{\pi }_{\max }}^{{\pi }_{\max }}}\mathrm{d}{s}_{||}{\xi }^{\mathrm{S}}(s_{||},s_{\perp })={\displaystyle {\int }_{-{\pi }_{\max }}^{{\pi }_{\max }}}\mathrm{d}{r}_{||}{\xi }^{\mathrm{R}}(r(r_{||},r_{\perp })),\end{array}$$(3)

where the equality holds if π_max is sufficiently large for correlations to be negligible beyond it. We use superscript S and R to denote redshift space and real space, respectively, and $r({\mathrm{r}}_{||},r_{\perp })\equiv \sqrt{{\mathrm{r}}_{||}^2+r_{\perp }^2}$. Importantly, the projected correlation function is independent of RSD effects under the distant-observer approximation, meaning we can model ξ^R without having to account for peculiar velocity, which can be directly confronted with measurements from observations.

3. Constructing mock galaxy catalogs

In this section we briefly describe the N-body simulations and the catalogs of HOD galaxies used in the analyses.

3.1. The ABACUSSUMMIT simulation suite

ABACUSSUMMIT (Maksimova et al. 2021) is a suite of large, high-accuracy cosmological N-body simulations. In this work, we mainly focused on the base simulations of the suite with a volume of 2 (h⁻¹ Gpc)³ containing 6912³ with a mass resolution of 2 ⋅ 10⁹ h⁻¹ M_⊙. The simulation suite contains 150 simulation boxes spanning 97 cosmological models. They assume a flat space-time, and the Hubble constant H₀ is calibrated to match the cosmic microwave background acoustic scale θ_*. The 8 varying parameters are

$$\begin{array}{c}\hfill \mathcal{C}=\{{\omega }_{\mathrm{cdm}},{\omega }_{\mathrm{b}},n_s,{\sigma }_8,w_0,w_a,\mathrm{d}{n}_s/\mathrm{d}lnk,N_{\mathrm{eff}}\},\end{array}$$(4)

where ω_cdm = Ω_cdmh² and ω_b = Ω_bh² are the physical cold dark matter and baryon densities, σ₈ is the amplitude of the matter power spectrum, w₀ is the dark energy equation of state today, w_a governing the time evolution of the dark energy equation of state according to the CPL (Chevallier & Polarski 2001; Linder 2003) parameterization,

$$\begin{array}{c}\hfill w=w_0+w_a\frac{z}{1+z},\end{array}$$(5)

where w₀ = −1 and w_a correspond to a cosmological constant, n_s is the spectral tilt, dn_s/dlnk is the running of the spectral tilt and N_eff is the effective number of neutrino species.

The base simulations considered in this paper are organized into the following subsets:

• c000: Fiducial cosmology matching Planck 2018 (Planck Collaboration VI 2020) particularly the mean estimates of the TT,TE,EE+lowE+lensing likelihood chains. This cosmology contains 25 independent realizations, referred to as the fiducial cosmology in this paper.

• c001-c004: Secondary cosmologies, including a low ω_cdm, thawing dark energy, high N_eff and low σ₈.

• c100-c126: Linear derivative grid with simulation pairs with small positive and negative steps in eight-dimensional parameter space.

• c130-c181: Broad emulator grid, covering a wider range of the full cosmological parameter space. All simulations in this grid have the same phase seed.

Table 1 shows the fiducial value and the prior ranges of the cosmological parameters covered by the simulations.

The c001-c004 simulations contain six independent realizations each, but we limited ourselves to a single realization to reduce data storage requirements. In addition to the base simulations, we used the ABACUSSMALL simulations to estimate the covariance matrix used for parameter inference analysis in Sect. 5. These contain 1687 independent realizations of the fiducial cosmology with 1728³ particles in a (0.5 h⁻¹ Gpc)³ volume box.

The group finding in each simulation is done on the fly, using a hybrid friends-of-friends spherical overdensity (SO) algorithm named COMPASO (Hadzhiyska et al. 2021). This paper only considers halo catalogs from snapshots at a single redshift of z = 0.25. As described in Hadzhiyska et al. (2021), there are three types of halos in the simulation suite, and we limited our analysis to the main halo catalogs, named L1. The resulting catalogs are then populated by galaxies, following the procedure described in the following subsection.

3.2. Modeling the galaxy-halo connection

To construct galaxy catalogs from the ABACUS halos, we adopted the HOD model introduced by Zheng et al. (2005), where the number of galaxies in a halo is determined by the halo mass alone, meaning that we ignored galaxy assembly bias. We parametrized central galaxies according to a Bernoulli distribution and satellites according to a Poisson distribution, using functional forms that were developed to match the clustering of luminous red galaxies.

The mean occupation number of central galaxies is parametrized as

$$\begin{array}{c}\hfill ⟨N_c⟩(M|\mathcal{G})=\frac{1}{2}(1+\mathrm{erf}\left(\frac{logM-log{M}_{\min }}{{\sigma }_{logM}}\right)),\end{array}$$(6)

where erf(x) is the error function, M_min represents the typical halo mass required to host a central galaxy, and σ_log M determines the width of the distribution.

The mean occupation number of satellites is given as

$$\begin{array}{c}\hfill ⟨N_s⟩(M|\mathcal{G})=⟨N_c⟩(M|\mathcal{G}){\left(\frac{M-\kappa M_{\min }}{M_1}\right)}^{\alpha },\end{array}$$(7)

where we make the common assumption that satellite galaxies only reside in halos that already contain a central galaxy. Here, κM_min is the minimum halo mass required to host a satellite galaxy, M₁ can be interpreted as the typical halo mass required to host a satellite galaxy, and α is a power-law index governing the satellite occupation’s mass dependence.

Our HOD model comprises the five HOD parameters. However, we can fix one parameter via the mean galaxy number density,

$$\begin{array}{c}\hfill {\overline{n}}_g={\displaystyle \int \mathrm{d}}M\frac{\mathrm{d}n}{\mathrm{d}M}[⟨N_c⟩(M|\mathcal{G})+⟨N_s⟩(M|\mathcal{G})],\end{array}$$(8)

where $\frac{\mathrm{d}n}{\mathrm{d}M}$ is the halo mass function. Given a set of 𝒢, we can replace one HOD parameter with ${\overline{n}}_g$ by inverting Eq. (8). Instead of M_min, we used ${\overline{n}}_g$ as a free parameter in our model since we have strong prior information about it from observation data. Hence, the free model parameters considered in this work is

$$\begin{array}{c}\hfill \mathcal{G}=\{{\overline{n}}_g,{\sigma }_{logM},M_1,\kappa ,\alpha \}.\end{array}$$(9)

We emphasize that this HOD model is highly simplistic, and several extensions have been developed over the years to accurately incorporate more complex properties of galaxy-halo connection supported by modern observations (e.g., Yuan et al. 2022). Nonetheless, the base model of Zheng et al. (2005) in HOD-related analyses remains widely used and is well tested, thus serving as a natural starting point to construct an emulator-based model for predicting galaxy clustering. Importantly, two-point statistics have intrinsic information loss, and we removed RSD effects by considering the projected correlation function. The base HOD model provides a benchmark for subsequent studies with more complex HOD parameterizations when additional summary statistics are included.

We followed Cuesta-Lazaro et al. (2023) in constructing mock galaxy catalogs, adopting R_200m as our halo radius definition, which is the radius of a sphere where the average enclosed density is 200 times the mean matter density of the Universe ρ_m0. The halo mass is given by the number of particles within the halo and is used to compute the halo radius. We note that our halo definitions deviate from those of COMPASO. In future analyses, we intend to fully integrate the COMPASO halo definitions to construct galaxy catalogs. Comparing results obtained from the two halo definitions will serve as a valuable test to determine the model independence of the results presented in this work, thereby ascertaining whether the underlying cosmology is indeed being probed.

Central galaxies are located at the host halos’ center, and satellite galaxies are distributed according to an Navarro-Frenk-White profile, u_NFW(r|c(M)) (Navarro et al. 1997). For the halo concentration, c(M), we used the median concentration-mass relation from Diemer & Joyce (2019). We assumed the angular distribution of the satellite galaxies to be isotropic.

To model RSD effects, we added velocity dispersions to the populated galaxies. By assumption, central galaxies follow their host halo and are assigned their velocity. Satellite galaxies are assigned the same velocity as the host halo, with an additional offset randomly drawn from a Gaussian distribution with mean zero and a standard deviation determined by the halo’s velocity dispersion. From the velocity dispersion of galaxies we get their redshift space position according to Eq. (2).

4. Emulating the galaxy two-point correlation function

To predict the projected correlation function given a set of cosmological parameters, 𝒞 and HOD parameters, 𝒢, we constructed a neural network emulator to predict ξ^R(r, 𝒞, 𝒢), from which w_p is obtained via Eq. (3). Emulators are computational models that efficiently approximate the behavior of a system, essentially acting as high-dimensional interpolators. In this section we describe the general setup of the neural network and its main prerequisites before introducing how mock galaxy catalogs are constructed from which the observable is computed.

To emulate the 2PCF, we employed fully connected neural networks, which approximate a function f such that

$$\begin{array}{c}\hfill \mathit{y}=f(\mathit{x}|\mathit{\theta }),\end{array}$$(10)

where x represents the features of the dataset, y the target values and θ the trainable parameters of the network.

The activation function of the network is the Gaussian error linear unit (GELU; Hendrycks & Gimpel 2016), defined as

$$\begin{array}{c}\hfill \mathrm{GELU}(x)=x·\mathrm{\Phi }(x)=x(1+\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)),\end{array}$$(11)

where x is the output of the network’s previous layer, and Φ(x) is the cumulative distribution function of the standard normal distribution. The GELU activation function is a smooth approximation of the Rectified Linear Unit (ReLU) activation function (Agarap 2019). GELU exhibits more curvature and nonlinearity, which may allow it to better approximate complex functions. Another advantage is that it avoids the vanishing gradient problem.

The optimal f is defined by the set of parameters θ, which minimizes a loss function. In this work, we adapted the L₁ norm

$$\begin{array}{c}\hfill L_1(\mathit{\theta })=\frac{1}{N}{\displaystyle \sum _i^N}|y_i^{\mathrm{pred}}(\mathit{\theta })-y_i^{\mathrm{true}}|,\end{array}$$(12)

where y_i^true is the true value of the observable and y_i^pred(θ) is the network prediction. We updated the network parameters using the Adam optimizer (Kingma & Ba 2017).

For the learning rate, we used a value of 0.01 initially and used a learning rate scheduler to reduce the learning rate by a factor of 10 if the validation loss did not improve after 20 epochs until a minimum value of 10⁻⁶ was reached. This simple scheme allows for increased training efficiency in the beginning when gradients are steep, followed by increasingly finer tuning of the parameters as the loss landscape flattens.

4.1. Training the neural network

We trained the neural network to perform the mapping:

$$\begin{array}{c}\hfill {log}_{10}(\xi (r))=f(\mathcal{C},\mathcal{G},r),\end{array}$$(13)

where 𝒞 and 𝒢 are the input parameters given in Eqs. (4) and (9), respectively, and r is the average separation distance between galaxy pairs in a given bin. Additionally, we applied standard normal scaling of the targets and features according to the transformation

$$\begin{array}{c}\hfill {\mathit{x}}_i\to \frac{{\mathit{x}}_i-⟨{\mathit{x}}_i⟩}{\sqrt{\mathrm{Var}({\mathit{x}}_i)}}·\end{array}$$(14)

To facilitate the generalization of the network, we employed weight decay and dropout. The former reduces the weights after each iteration, preventing significant emphasis from being placed on the noise of outliers, essentially working by subtracting an L₂ penalization term λ||θ||₂ from Eq. (12). Dropout counteracts co-adaption, in which each network neuron depends on specific activations of other neurons, by randomly setting neurons to 0 with a probability p. Given that we have a restricted number of samples with varying cosmological parameters, dropout is an effective method for enhancing generalization (Srivastava et al. 2014). To further prevent overfitting, the training was stopped if the validation loss did not improve after 50 epochs or when a maximum of 500 epochs was reached.

4.2. Generating mock-galaxy catalogs

Our choice of fiducial HOD parameters is principally based on the analysis of Cuesta-Lazaro et al. (2023) who applied the same HOD model on the DARKQUEST simulations (Nishimichi et al. 2019) where the minimum halo mass is 10¹² h⁻¹ M_⊙. In the context of these HOD parameters, ABACUS halos of mass M < 10¹² h⁻¹ M_⊙ are improbable to host galaxies according to the occupation distributions in Eqs. (6) and (7). This consideration, combined with the benefits of reduced computational cost and data storage requirements of the generated mock-galaxy catalogs, prompts us to restrict our analysis to halos of mass M ≥ 10¹² h⁻¹ M_⊙. There are substantial low-mass halos in the simulation suite, which means they likely contribute with additional galaxies despite have a low probabilities of hosting any. Our approach is overly simplistic, and an alternative approach to reducing the number of halo samples is to randomly down-sample halos according to a probability distribution that depends on the halo mass, as considered by Yuan et al. (2022, Eq. 13).

Instead of M_min, we used ${\overline{n}}_g$ as a free parameter in our model since we have strong prior information about it. For the remaining HOD parameters, we adopted the fiducial values used by Cuesta-Lazaro et al. (2023, Table 3) and chose variations of ±10% as our prior ranges. Our HOD parameter’s fiducial values and prior ranges are shown in Table 2.

We randomly assigned 80% of the c100-c126 and c130-c180 cosmologies to the training sets and the remaining 20% to the validation sets. The c000-c004 were assigned to the test sets. Degeneracies between cosmological and HOD parameters are expected, and to break these degeneracies, we constructed numerous catalogs for each simulation. Breaking these degeneracies requires using numerous HOD parameter combinations for each cosmology. We therefore generated 500 sets of HOD parameters for the training set and 100 parameter sets each for the test and validation sets. To evenly generate samples of HOD parameters within the prior ranges in Table 2, we generated five-dimensional Latin hypercubes using the Python module SMT (Saves et al. 2024). The cosmological parameters space covered by the ABACUSBASE simulations are shown in Fig. 1.

Fig. 1. Cosmological parameters used in the training, testing, and validation datasets. Note that the fiducial parameter values are included in the test dataset. The large number of fiducial values for some parameters is due to the linear derivative varying two parameters at a time and the fixed parameter subsets of the broad emulator grid.

As mentioned in Sect. 3.2 we estimated the appropriate value of M_min by inverting Eq. (8) to generate galaxy catalogs with galaxy number densities matching the desired value. From the generated galaxy catalogs, we stored the coordinates of all galaxies in real space, (x,  y,  z), and the redshift space coordinates s_z as we isolated RSD effects to the z-axis.

4.3. Computing the 2PCF and its projection

We measured the real space 2PCF with the Python module PYCORR¹, utilizing the CORRFUNC engine (Sinha & Garrison 2020; Majumdar & Arora 2019). Measurements were made using 40 logarithmically spaced separation bin edges of r ∈ [0.01,  5) h⁻¹ Mpc and 75 linearly spaced bin edges of r ∈ [5,  150] h⁻¹ Mpc. These bin edges are chosen to cover a range of scales, ensuring a relatively high resolution on nonlinear scales where we expected significant clustering and lower resolution on larger scales where galaxy pair counts are expected to be sparse. We retrieved the average separation distance between galaxy pairs in each bin from the measurements, which we used as input for the emulator.

Initial emulation tests revealed substantial errors in the predicted clustering on the largest and smallest scales, despite excluding samples with ξ ≤ 10⁻⁷, as shot noise was prominent. To mitigate this, we opted to restrict the analysis to samples within the range r ∈ (0.1,  105) h⁻¹ Mpc, containing only samples with ξ^R(r) > 10⁻⁷.

To compute w_p(r_⊥) from ξ^R(r), we used 40 logarithmically spaced bin edges of r_⊥ ∈ [0.5,  40] h⁻¹ Mpc. Unless otherwise stated, we use this binning for all computations of w_p throughout this paper. We then created a spline of ξ^R(r), which we integrated along r_|| ∈ [0, 105] h⁻¹ Mpc according to Eq. (3), setting ξ^R = 0 for values where r(r_⊥, r_||) extended beyond the training range considered.

5. Recovering cosmological parameters

In this section we demonstrate our emulator model’s ability to predict cosmological parameters from a mock signal of the projected correlation function. In the future, we aim to apply our model to observational data from the ongoing Euclid survey (Laureijs et al. 2011). In this work, we assessed the model’s ability to recover cosmological parameters using mock data, which provides a controlled environment for evaluating its limitations and constraining power, setting a benchmark against unknown data. Importantly, the mock data help us identify other summary statistics to incorporate into our model before applying it to real data.

We performed Markov chain Monte Carlo (MCMC) parameter inference to sample the posterior distributions in our 13-dimensional parameter space. We began by defining the posterior distribution. Up to a constant, it is given as

$$\begin{array}{c}\hfill p(\theta |\mathcal{D})\propto \mathcal{L}(\mathcal{D}|\theta )p(\theta ),\end{array}$$(15)

where θ = {𝒞,𝒢} are the parameters we want to estimate, 𝒟 is the mock data, p(θ | 𝒟) is the posterior of the parameters given the data and p(θ) is the prior distribution of the parameters.

The likelihood function ℒ(𝒟 | 𝒞, 𝒢) describes how well our model, specified by our parameters, explains the observed data. We assumed that the likelihood follows a multivariate Gaussian distribution, as is common when retrieving cosmological parameters from galaxy clustering observables (Cuesta-Lazaro et al. 2024; Miyatake et al. 2022), and compute the log-likelihood as

$$\begin{array}{cc}\hfill ln\mathcal{L}(\mathcal{D}|\theta )=& -\frac{1}{2}{\displaystyle \sum _{r_{\perp }^i,r_{\perp }^j}}[w_p^{\mathrm{data}}(r_{\perp }^i)-w_p(r_{\perp }^i|\theta )]\hfill \\ \hfill & \times {\mathit{C}}^{-1}(w_p^{\mathrm{data}}(r_{\perp }^i),w_p^{\mathrm{data}}(r_{\perp }^j))\hfill \\ \hfill & \times [w_p^{\mathrm{data}}(r_{\perp }^j)-w_p(r_{\perp }^j|\theta )]\hfill \\ \hfill & +\mathrm{const}.,\hfill \end{array}$$(16)

where w_p^data(r_⊥) is the mock data vector, w_p(r_⊥ | θ) is the emulator’s prediction given the model parameters and C⁻¹ is the inverse of the covariance matrix.

We used the 25 fiducial c000 simulation boxes to construct the mock data vector, using the fiducial HOD parameters listed in Table 2. We used CORRFUNC (Sinha & Garrison 2020) to estimate the projected correlation function in redshift space. The mock data signal was taken as the average w_p(r_⊥) from the 25 samples, shown in Fig. 2.

Fig. 2. Mock data vector wpdata(r⊥) (red) computed from averaging over the individual wp from the 25 fiducial cosmologies (blue). The error bars show the standard deviation of the individual wp samples.

From the 25 data vector catalogs we also computed ξ^R over the emulator’s full separation training range, r ∈ (0.1,  105) h⁻¹ Mpc. The average separation vector from the 25 samples as our fixed separation vector to predict ξ^R(r | θ) during parameter inference, from which w_p(r_⊥ | θ) is obtained.

We estimated the covariance matrix using the 1687 ABACUSSMALL simulation boxes, generating one mock-galaxy catalog for each box using the fiducial HOD parameters in Table 2. A cutoff mass is not considered for these simulations, as each box contains significantly fewer halos than the ABACUSBASE boxes. The projected correlation function for each catalog is computed according to the same prescription used to compute the mock data vector. The covariance matrix is scaled by a factor of 1/64 to match the volume of the base simulations. Figure 3 shows the normalized covariance matrix. The covariance is strictly positive, and off-diagonal terms are dominated by signals with a high signal-to-noise ratio, resulting from the large number of samples considered. The small-scale bins are primarily independent, as they are dominated by shot noise of the galaxy catalogs, whereas the large-scale bins are correlated due to cosmic variance.

Fig. 3. Normalized covariance matrix for wp(r⊥), computed from the 1687 fiducial cosmology simulations in the ABACUSSMALL dataset.

Sampling is performed using the EMCEE python package (Foreman-Mackey et al. 2013), employing the differential evolution algorithm (Braak 2006). The MCMC chain is assigned initial parameter values of 𝒞_fid and 𝒢_fid, where a small Gaussian perturbation of 𝒩(0, 10⁻³) is added to each parameter. Samples are drawn from a flat prior distribution of 𝒞 and 𝒢, with prior ranges given in Tables 1 and 2, respectively.

In addition to the inference analysis, we sampled the posterior distribution of the cosmological parameters with HOD parameters fixed to their fiducial value and vice versa, which helped us identify and address degeneracies between the two parameter sets. This analysis is vital to identify potential supplementary statistics to add to our model to better constrain HOD parameters. Such tailored analysis is unsuitable for real data due to underestimation and bias in error measurements, underscoring the mock dataset’s importance during model development.

6. Results

6.1. Emulator performance

The emulator’s prediction of ξ^R is shown in Fig. 4, with one random sample from every cosmology in the test dataset shown in the upper panel. In the bottom panel, we show relative errors of the individual samples from the upper panel and include the mean and median errors computed from the test dataset. The emulator accuracy increases as the extreme separation distances are approached, where Poisson noise becomes increasingly dominant. On intermediate scales, the emulator achieves sub-percent accuracy with relative errors below 1% for 1 ≲ r/(h⁻¹ Mpc)≲30.

Fig. 4. Upper panel: Example set of ξR sampled from 29 cosmologies in the test and validation datasets, comparing data measurements with emulator predictions. We scale ξR with r to emphasize the accuracy on large scales. Lower panel: Relative error of the prediction samples from the upper panel (gray lines) and the mean and median relative error from all samples in the test data.

Figure 5 shows the emulator prediction of w_p in the upper panel, using the same test samples as in Fig. 4. The data points correspond to direct measurements of w_p using the redshift space coordinates of galaxies. Relative errors from each sample are shown in the lower panel, including mean and median errors computed from all test simulations. Compared to the errors of ξ^R, the emulator predicts w_p more accurately, with errors below 2% for r_⊥ ≲ 20 h⁻¹ Mpc. The projection results in an overall smoothing of the relative errors of w_p compared to ξ^R, with reduced errors toward the extreme ends and increased errors on intermediate scales. This indicates that some information is lost during the projection of ξ^R.

Fig. 5. Upper panel: Example set of wp sampled from 29 cosmologies in the test and validation datasets, with solid lines corresponding to the integrated samples of ξR in Fig. 4. Data points show direct measurements of wp in redshift space. Lower panel: Relative difference between wp from direct measurements and from integrating ${\xi }_{\mathrm{pred}}^{\mathrm{R}}$.

The largest emulator errors corresponded to specific HOD parameter combinations that produced noisy samples of ξ^R on large scales where ξ^R dropped by orders of magnitude from one measurement to the next. These outliers occurred after limiting the separation scales to r ≤ 105 h⁻¹ Mpc. Combined with the prominent shot noise in ξ^R at r ≳ 60 h⁻¹ Mpc, this suggests that the bin resolution was too narrow on large scales.

We note that the c002 simulation has the highest value of w₀ of all cosmologies considered, and is thus outside the training range of the emulator. Nonetheless, relative error from this simulation was lower than the mean and median error from test simulations. Additionally, we did not ensure that the extreme values of the generated HOD parameters belonged to the training set, which caused some to be found in the testing and validation set.

6.2. Parameter inference

Figure 6 shows the posterior distributions of the cosmological parameters, marginalized over the HOD parameters. The posterior of baryonic density ω_b is poorly constrained using the projected 2PCF data alone, and therefore not shown in the plot. This is also a consequence of using an overly simplistic HOD model and the absence of the information of baryonic acoustic oscillations. Several parameters have wide marginalized distributions, indicating that our model cannot put significant constraints on them from w_p at z = 0.25 over the separation scale considered.

Fig. 6. Recovery test of the fiducial cosmology after marginalizing over the HOD parameters (in blue) and keeping the HOD parameters fixed (in red).

The posterior distributions of ω_cdm and σ₈ suggest a tighter constraint on these parameters since the two parameters directly impact structure formation. There appears to be no considerable degeneracy between ω_cdm and the HOD parameters. On the other hand, σ₈ shows distinct levels of degeneracy with the HOD parameters since its effect on the 2PCF can easily be compensated by changes in HOD parameters. Similar results were obtained by Miyatake et al. (2022, Fig. 9), who were able to constrain Ω_m and σ₈ well using w_p only using the DARKQUEST data.

We note that one of the equation of state parameters, w_a, is poorly constrained by our model. This can be explained by the fact that our analysis is restricted to a single redshift, we cannot expect our model to appropriately account for dynamic dark energy. However, Cuesta-Lazaro et al. (2024, Fig. 6) obtains similar constraints for these parameters using ξ^S, but their constraints are significantly improved when including signals from density-split clustering, suggesting that the two parameters could be constrained with a single redshift only.

The HOD parameter posterior distributions are shown in Fig. 7. The HOD parameters that are well constrained are ${\overline{n}}_g$ and M₁, which are the two parameters most strongly affecting the central and satellite galaxy distributions. Both parameters show some degeneracy with the cosmological parameters, most notably in the marginalized distributions of ${\overline{n}}_g$. This is expected, as cosmological parameters affect the resulting ${\overline{n}}_g$ through the halo mass function. The parameters σ_{log M} and α are relatively poorly constrained. Their posteriors, when varying HOD parameters only, suggest that wider prior ranges should be considered, as in Cuesta-Lazaro et al. (2024, Table 1). Our results are consistent with previous studies, such as Zentner et al. (2019, Fig. 1), in which HOD parameters are poorly constrained by w_p alone.

Fig. 7. Marginalized posterior of the HOD parameters. The true parameters in the recovery test are shown in black.

To further validate our pipeline, we performed additional recovery tests using two different cosmologies. The results are presented in Fig. 8. We find that the true cosmological parameters are consistently recovered within 2σ for both cosmological models. This demonstrates the robustness of our method across different cosmologies and supports the reliability of our parameter constraining procedure.

Fig. 8. Recovery tests for two different fiducial cosmologies selected from the test dataset. The true parameter values are indicated by the colored horizontal and vertical dashed lines.

7. Discussion

Our analysis of parameter degeneracies is mainly limited to degeneracies between the cosmological and HOD parameters. To accurately discern the impact on a specific parameter, it is essential to consider the intrinsic constraining power that the parameter in question exerts on the observable. Our analysis of these degeneracies is hindered by the model’s limited ability to tightly constrain several parameters. As a result, further extensions to the model are necessary to adequately discriminate the effects of the two parameter sets. Moreover, degeneracies within each parameter set are not addressed in our current study.

A significant limitation of our recovery test is the primitive estimation of the covariance matrix. Most notably, emulator errors are not accounted for. As illustrated in Fig. 5, the predictions of our emulator are less accurate at extreme scales, which is a factor that is not incorporated into the likelihood model in Eq. (16). Furthermore, the errors among individual measurements of w_p used to construct the mock data vector, as shown in Fig. 2, are neglected. Developing a more realistic covariance matrix is the essential first step toward applying this emulator-based approach to real data. Consequently, we are unable to draw definitive conclusions about the true constraining power of our model. After refining the current likelihood estimation of our model, a second clustering statistic should be added before further explorations of the information content of our inference analysis is considered, such as the scale dependence of parameter’s constraints.

Bias parameters are among the critical ingredients missing in the base model of Zheng et al. (2005). Both halo assembly bias, referring to halos of equal mass exhibiting different clustering properties, and galaxy assembly bias, referring to variations in the number of galaxies hosted by halos of identical mass, have been omitted in our analysis. Cuesta-Lazaro et al. (2023, Fig. 10), for instance, shows how including environmental-based assembly bias can induce a shift in the constraint of σ₈ more than 1σ from their fiducial value. Other effects, such as velocity biases, are often introduced to allow a galaxy’s peculiar velocity to be slightly offset from that of its host halo, which is supported by observational evidence (Yuan et al. 2022). In future analyses, we plan to construct refined galaxy catalogs using the highly efficient PYTHON package ABACUSHOD (Yuan et al. 2022), which includes a large variety of modern HOD variations and is specifically designed to work with the ABACUSSUMMIT simulation suite.

Several works in the literature also emulate galaxy clustering statistics using different simulation suites, such as the AEMULUS project (Zhai et al. 2019, 2023; Storey-Fisher et al. 2024) and the Mira-Titan project (Kwan et al. 2023). These simulations cover a similar range of cosmological parameter space as the ABACUS simulations, but feature different extensions to the baseline Λ cold dark matter model, such as massive neutrinos and dynamical dark energy. In contrast to previous works, which typically use Gaussian process regression, we adopted fully connected neural networks for our emulator construction. Neural networks are particularly well suited for handling large datasets, unlike Gaussian process regression, which scales as 𝒪(N³), with N being the size of the training dataset.

8. Conclusions

Emulators provide an efficient method for bypassing the prohibitively expensive computational demands of running full N-body simulations to model galaxy clustering statistics on nonlinear scales. We have introduced a model designed to predict the real-space correlation function of galaxies over the range r ∈ [0.1,  105] h⁻¹ Mpc. The key strength of this approach lies in its computational efficiency, which enables a near-instantaneous estimation of the correlation function given a set of cosmological parameters. This was achieved by constructing galaxy catalogs using a straightforward HOD model applied to halos from the ABACUSSUMMIT simulation suite.

Our emulator predicts the projected correlation function over a wide range of scales, with sub-percent accuracies on scales 1 ≲ r_⊥/(h⁻¹ Mpc)≲10. Using this model, we conducted inference analyses to constrain key cosmological parameters, particularly ω_cdm and σ₈. While these analyses demonstrate the model’s capability, the limited information provided by the projected correlation function alone prevents us from placing stronger constraints on other cosmological parameters. This issue can be addressed by incorporating the information from redshift-space clustering, which we plan to explore in future work. Furthermore, our analysis reveals that the simplified covariance matrix employed does not account for emulator errors, highlighting the need for further refinement of our methodology before drawing definitive conclusions. Despite these limitations, our results are consistent with previous studies regarding the constraining power of the projected correlation function.

Future work in developing an emulator-based model for large-scale structure analysis will focus on incorporating additional summary statistics to enhance the model’s ability to probe underlying cosmological information, such as higher-order clustering statistics, density-split clustering (Cuesta-Lazaro et al. 2024; Paillas et al. 2024), and the marked correlation function (Storey-Fisher et al. 2024). Additionally, by employing a more sophisticated HOD model, we aim to generate galaxy catalogs that more accurately reflect the complexities of galaxy clustering. At present, our model cannot robustly marginalize over the HOD parameters, a limitation we intend to address. Once this is achieved, we plan to apply our model to observational data from the ongoing Euclid survey to further refine cosmological parameter estimates.

Acknowledgments

We thank the Research Council of Norway for their support and the resources provided by UNINETT Sigma2 – the National Infrastructure for High-Performance Computing and Data Storage in Norway. We are grateful to the editor and the anonymous reviewer for their valuable comments and suggestions.

References

All Tables

All Figures

	Fig. 1. Cosmological parameters used in the training, testing, and validation datasets. Note that the fiducial parameter values are included in the test dataset. The large number of fiducial values for some parameters is due to the linear derivative varying two parameters at a time and the fixed parameter subsets of the broad emulator grid.
In the text

	Fig. 2. Mock data vector wpdata(r⊥) (red) computed from averaging over the individual wp from the 25 fiducial cosmologies (blue). The error bars show the standard deviation of the individual wp samples.
In the text

	Fig. 3. Normalized covariance matrix for wp(r⊥), computed from the 1687 fiducial cosmology simulations in the ABACUSSMALL dataset.
In the text

	Fig. 4. Upper panel: Example set of ξR sampled from 29 cosmologies in the test and validation datasets, comparing data measurements with emulator predictions. We scale ξR with r to emphasize the accuracy on large scales. Lower panel: Relative error of the prediction samples from the upper panel (gray lines) and the mean and median relative error from all samples in the test data.
In the text

	Fig. 5. Upper panel: Example set of wp sampled from 29 cosmologies in the test and validation datasets, with solid lines corresponding to the integrated samples of ξR in Fig. 4. Data points show direct measurements of wp in redshift space. Lower panel: Relative difference between wp from direct measurements and from integrating ${\xi }_{\mathrm{pred}}^{\mathrm{R}}$.
In the text

	Fig. 6. Recovery test of the fiducial cosmology after marginalizing over the HOD parameters (in blue) and keeping the HOD parameters fixed (in red).
In the text

	Fig. 7. Marginalized posterior of the HOD parameters. The true parameters in the recovery test are shown in black.
In the text

	Fig. 8. Recovery tests for two different fiducial cosmologies selected from the test dataset. The true parameter values are indicated by the colored horizontal and vertical dashed lines.
In the text