Full Shape Cosmology Analysis from BOSS in configuration space using Neural Network Acceleration

Over the last decade, the study of the clustering of galaxies through Large Scale Structure (LSS) surveys has emerged as a crucial probe within precision Cosmology. Spectroscopic surveys, such as the Sloan Digital Sky Survey²²2www.sdss.org (SDSS) and the Dark Energy Spectroscopic Instrument³³3www.desi.lbl.gov (DESI) provide three-dimensional maps of the Universe, where angular positions and redshifts of millions of galaxies are measured with high accuracy. These maps constitute appropriate data sets for quantifying the clustering characteristics of galaxies, including correlation functions and other summary statistics, which then can be compared to models’ theoretical predictions.

The study of clustering statistics has primarily relied on two significant sources of cosmological information. First, Baryon Acoustic Oscillations (BAO) in the early Universe freeze up at the drag epoch, whose signature can be observed at later times in the correlation function as a well-distinguished peak around a scale of 150 Mpc. Second, peculiar velocities of galaxies contributes to the redshift that we measure adding to the Hubble flow component. Since this contribution occurs only along the line-of-sight direction, we observe an apparent anisotropic distortion in the matter distribution, and hence, galaxy statistics which otherwise would be isotropic become dependent on the angle of observation. This effect is known as Redshift Space Distortions (RSD).

RSD and BAO effects most of the relevant information about the correlation function of galaxies. Consequently, the SDSS-III BOSS [1, 2, 3] and SDSS-IV eBOSS [4, 5, 6] collaborations have chosen compressed methodologies for their standard analysis. In such approaches, the cosmological parameters of the matter power spectrum are fixed to fiducial values, and a set of parameters characterizing the BAO and RSD effects are explored. Given that RSD depends on the average velocity of galaxies, it is sensitive to the growth rate of structure, so one of the chosen parameters is $f\sigma_{8}$. Two additional degrees of freedom should be included to account for the distortions in the position of the BAO along and across the line-of-sight, which arise from potential mismatches between the fiducial and true cosmologies, that is, the Alcock-Paczyński effect [7].

On the other hand, the Effective Field Theory of LSS (hereafter simply EFT) [8, 9, 10, 11, 12, 13] built on top of Perturbation Theory [14] has been developed during the past years, and by now is currently used in analyses that confront theoretical models of the galaxy distribution directly to the data gathered by our telescopes. These methods are commonly known as full-modeling or full-shape analyses, and operate in a similar fashion that it has been done for the CMB over the years. Nowadays, these full-shape templates are used routinely to constrain cosmological parameters [15, 16, 17, 18, 19, 20], including higher-order statistics [21, 22] and even beyond $\Lambda$CDM models [23, 24, 25, 26, 27].

One of the primary advantages of the compressed methodology is its agnostic nature: the parameters it explores are relatively model-independent when compared with those obtained from a direct, full-shape analysis. However, generating full-shape theoretical templates of the power spectrum comes with significant computational costs, which, until recently, hindered our ability to perform cosmological parameter estimation. This is one of the reasons why the compressed methodology has been favored by some part of the community. But even if the full-shape analysis is expensive and model-dependent, it is capable of extracting more cosmological information from the power spectrum. Therefore, both the compressed and full-shape methodologies have their own merits, and there should be an incentive to pursue both approaches in parallel.

As we transition into a new era of cosmological surveys, the costs of full-shape models are set to increase even further. This is due to the unprecedented precision achieved in measurements on small scales of the correlation function (smaller than 50 $h^{-1}\,\textrm{Mpc}$), where the nonlinearities of perturbations exert a none negligible influence on halo distributions. Furthermore, at these small scales, the relationship between the clustering of observable astrophysical tracers and the underlying dark matter halos becomes complex. As a consequence, building accurate templates at these small scales might require the evaluation of even more complicated models, thereby introducing more intricate calculations and increasing the computational cost of an individual template. Finally, the analytical modeling of higher order statistics will pose new computational time challenges.

In recent years, a number of avenues have been developed to tackle these difficulties and going beyond the standard three parameters compress methodology and into more intricate models at this smaller scales. One possible road is to expand the compress approach by introducing a small subset of new free parameters that encompass most of the remaining relevant information in the power spectrum [28]. Another possibility is to perform full shape analyses encompassing various cosmological parameters. To achieve this, several optimizations in the computational methods used to construct theoretical templates have been developed, significantly reducing the computational cost of full shape analyses.

As a result, many groups have reanalyzed the BOSS and eBOSS data using full-shape modeling. The optimized methodologies employed for these analyses can be categorized into two groups: efficient theoretical templates of the power spectrum [e.g. 15, 16, 17, 29, 30, 31] or emulator techniques that learn how to reproduce expensive models [e.g. 32]. The consensus from all of these reanalysis methods is that the constraints on cosmological parameters like the Hubble constant $H_{0}$ are significantly tighter when using these improved approaches [33].

In configuration space, the number of analyses in the literature is reduced, mainly because the consensus is that direct-fits are more constrictive in Fourier space. Nevertheless, working directly in configuration space has its benefits, particularly when dealing with the well localized BAO peak, which in Fourier space becomes distributed across a wide range wave numbers. Adding to that, some of the observational systematic have different effects in Fourier and configuration space. Consequently there is an incentive to study both spaces simultaneously. In ref.[30] the main analysis is performed in Fourier space but the authors have also worked out the correlation function as a consistency check. On the other hand, the work of [19] is devoted to perform fitting only to the correlation function using the PyBird code. The main difference with our approach is that we work from the beginning in a Lagrangian framework, and get directly the correlation function without the necessity of obtain first the power spectrum (plus adding infrared resummations) as an intermediate steps, to at the end Fourier Transform the results to the obtain the correlation function.

The number of free parameters to explore in full shape analysis is generally large when compared to the standard approach. With this in mind, considerable efforts have been made in building efficient sampling methods [34, 35, 36] which reduce the number of model evaluations required to run Markov Chain Monte Carlo (MCMC) explorations. However, the models will most likely continue to increase in complexity in the next years, since higher loop contribution or higher order statistics can be considered in the analyses. Therefore, there is an incentive to reduce the computational cost of generating these theoretical templates.

Machine Learning algorithms like Neural Network have been successfully used to drastically reduce the evaluation times of complex models [e.g. 37]. These techniques use datasets of pre-computed templates at different points within the parameter space and learn how to reproduce them. When trained correctly, neural networks can reproduce fairly complex models with an error smaller than the precision needed by LSS surveys. Also, given that neural networks are not local interpolators, in principle the errors in their predictions are not as strongly dependent on the distance of the nearest point within the training set, as methodologies like Gaussian processes emulators would be.

In this work, we model the redshift space correlation function up to one-loop perturbation theory using a Gaussian Streaming model [38, 39, 40, 41, 12] in combination with Effective Field theory (EFT). Throughout this work, we refer to our modeling as EFT-GSM. To implement our model, we release a code⁴⁴4https://github.com/alejandroaviles/gsm that uses a brute force approach, but still can compute the correlation function in $\mathcal{O}(1\,\text{sec})$ time, as described in section 3.1. However, the number of evaluations required to build a convergent MCMC chain for our baseline analysis is large, and this process can take a considerably amount of time. Therefore, we also built a neural network emulator to accelerate the computation of individual templates, reducing the running time of an MCMC chain from a few tens of hours to around 60 minutes (using the same computing settings) and below 20 minutes when we run the code in parallel depending on the cluster settings.

We utilize our methodology to reanalyze the BOSS DR12 LRG data obtaining the tightest constraints on the $\Lambda$CDM parameters using the 2-point correlation function alone. Throughout this study, we pursue two primary objectives. The first is to bring forward the potential of our full shape modelling approach in configuration space for extracting cosmological information when compared to its counterpart in Fourier space. The second objective is to demonstrate that neural network surrogate models can be used safely to optimize cosmological analysis, leading to significant savings in both time and computational resources without sacrificing accuracy when analyzing real data.

We have also tested extended regions in parameter space compared with the baseline analysis to include scenarios with prior configurations beyond Planck [42] and Big Bang Nucleosynthesis (BBN) [43] on parameters $\Omega_{b}$ and $n_{s}$ that will serve us to explore the potential of LSS observables standalone. This serves as a test of the full methodology in highly more degenerated scenarios than the baseline and with aim to prove the viability of the use of neural network in this larger and more complex parameter space.

This paper is organized as follows, we begin in section 2, introducing the data from the BOSS collaboration that we analyze here, as well as introducing a set of different mock simulations that we use for testing our methodology and building the covariance matrices required for our likelihood estimations. Then in section 3 we introduce the EFT-GSM model that we utilize to construct our theoretical templates. Then, in section 4 we introduce our neural network methodology, which is used as a surrogate model instead of the EFT-GSM model, we also quantify how much efficiency is gained with this. In section 5 we describe the fitting methodology, including a brief description of full shape fits. Here we also emphasise our parameter space and the priors that we impose on each parameter. Section 6 presents the validation of the methodology with high precision mocks and section 7 presents the results of our baseline analysis and how it compares with published alternative analysis of BOSS. We also include a subsection with results expanding the priors for cosmological parameters usually constrained independently by other observables and a subsection exploring the information content in the multipoles.

In this work we adopt a Lagrangian approach, on which we follow the trajectories of cold dark matter particles with initial position $\bm{q}$ through the map

$$\bm{x}(\bm{q},t)=\bm{q}+\mathbf{\Psi}(\bm{q},t),$$

(3.1)

where $\mathbf{\Psi}$ is the Lagrangian displacement field. The observed positions of objects are distorted by the Doppler effect induced by their peculiar velocities relative to the Hubble flow, $\bm{v}=a\dot{\bm{x}}=a\mathbf{\dot{\Psi}}$. That is, for a tracer located at a comoving real space position $\bm{x}$, its apparent redshift-space position becomes $\bm{s}=\bm{x}+\bm{u}$, with along the line-of-sight “velocity” $\bm{u}$

$$\bm{u}\equiv\hat{\bm{n}}\frac{\bm{v}\cdot\hat{\bm{n}}}{aH}=\hat{\bm{n}}\frac{\mathbf{\dot{\Psi}}\cdot\hat{\bm{n}}}{H},$$

(3.2)

where we are using the distant observer approximation on which the angular observed direction of individual galaxies $\hat{\bm{x}}_{i}$ are replaced by a single line-of-sight direction $\hat{\bm{n}}$, which is representative to the sample of observed objects. The map between Lagrangian coordinates and redshift-space Eulerian positions becomes

$$\bm{s}=\bm{q}+\mathbf{\Psi}+\hat{\bm{n}}\frac{\mathbf{\dot{\Psi}}\cdot\hat{\bm{n}}}{H}.$$

(3.3)

The correlation function of tracer counts in redshift space is given by the standard definition

$$1+\xi_{s}(\bm{s})=\langle\big{(}1+\delta_{s}(0)\big{)}\big{(}1+\delta_{s}(\bm{s})\big{)}\rangle$$

(3.4)

where the tracer fluctuation in redshift space is $\delta_{s}(\bm{s})$. Now, the conservation of number of tracers $X$, that reads $\big{[}1+\delta_{s}(\bm{s})\big{]}d^{3}s=\big{[}1+\delta_{X}(\bm{x})\big{]}d^{3}x$, yields the redshift-space correlation function [39]

$$1+\xi_{s}(\bm{s})=\int\frac{d^{3}k}{(2\pi)^{3}}d^{3}r\,e^{i\bm{k}\cdot(\bm{s}-\bm{r})}\Big{[}1+\mathcal{M}(\bm{k},\bm{r})\Big{]},$$

(3.5)

where $\bm{r}=\bm{x}_{2}-\bm{x}_{1}$ and $\bm{s}=\bm{s}_{2}-\bm{s}_{1}$. The density weighted pairwise velocity generating function is

$$1+\mathcal{M}(\bm{J},\bm{r})=\left\langle\big{(}1+\delta_{X}(\bm{x}_{1})\big{)}\big{(}1+\delta_{X}(\bm{x}_{2})\big{)}e^{-i\bm{J}\cdot\Delta\bm{u}}\right\rangle,$$

(3.6)

where $\Delta\bm{u}=\bm{u}(\bm{x}_{2})-\bm{u}(\bm{x}_{1})$.

The generating function is now expanded in cumulants $\mathcal{C}$. That is,

$$\mathcal{Z}(\bm{J},\bm{r})\equiv\log\Big{[}1+\mathcal{M}(\bm{J},\bm{r})\Big{]}=\sum_{n=0}^{\infty}\frac{(-i)^{n}}{n!}J_{i_{1}}\cdots J_{i_{n}}\mathcal{C}^{(n)}_{i_{1}\cdots i_{n}}(\bm{r})$$

(3.7)

where in the second equality we use the Taylor series of $\mathcal{Z}$ about $\bm{J}=0$. The cumulants are then obtained by

$$\mathcal{C}^{(n)}_{i_{1}\cdots i_{n}}(\bm{r})=i^{n}\frac{\partial^{n}\mathcal{Z}(\bm{J},\bm{r})}{\partial J_{i_{1}}\cdots\partial J_{i_{n}}}\Bigg{|}_{\bm{J}=0}.$$

(3.8)

Then,

$$1+\xi_{s}(\bm{s})=\int\frac{d^{3}k}{(2\pi)^{3}}d^{3}x\,e^{i\bm{k}\cdot(\bm{s}-\bm{r})}\exp\left[\sum_{n=0}^{\infty}\frac{(-i)^{n}}{n!}k_{i_{1}}\cdots k_{i_{n}}\mathcal{C}^{(n)}_{i_{1}\cdots i_{n}}(\bm{r})\right],$$

(3.9)

On the other hand, the generating function can be alternatively expanded in moments $\Xi^{(n)}$ as follows

	$\displaystyle\Xi^{(n)}_{i_{1}\cdots i_{n}}(\bm{r})$	$\displaystyle=i^{n}\frac{\partial^{n}\,}{\partial J_{i_{1}}\cdots J_{i_{n}}}\Big{(}1+\mathcal{M}(\bm{J},\bm{r})\Big{)}\Bigg{|}_{\bm{J}=0}$
		$\displaystyle=\left\langle\big{(}1+\delta_{X}(\bm{x}_{1})\big{)}\big{(}1+\delta_{X}(\bm{x}_{2})\big{)}\Delta u_{i_{1}}\cdots\Delta u_{i_{n}}\right\rangle,$		(3.10)

which lead us to relations between cumulants and moments

	$\displaystyle\mathcal{C}^{(0)}(\bm{r})$	$\displaystyle=\log[1+\xi(r)],$		(3.11)
	$\displaystyle\mathcal{C}^{(1)}_{i}(\bm{r})$	$\displaystyle=\frac{\Xi^{(1)}_{i}(\bm{r})}{1+\xi(r)}\equiv v^{\hat{n}}_{12,i},$		(3.12)
	$\displaystyle\mathcal{C}^{(2)}_{ij}(\bm{r})$	$\displaystyle=\frac{\Xi^{(2)}_{ij}(\bm{r})}{1+\xi(r)}-\mathcal{C}^{(1)}_{i}(\bm{r})\mathcal{C}^{(1)}_{j}(\bm{r})\equiv\hat{\sigma}^{2\hat{n}}_{12,ij}-v^{\hat{n}}_{12,i}v^{\hat{n}}_{12,j}=\sigma^{2\hat{n}}_{12,ij},$		(3.13)

where we introduced the pairwise velocity along the line of sight $v^{\hat{n}}_{12,i}$ and the pairwise velocity dispersion along the line of sight moment and cumulant, $\hat{\sigma}^{2\hat{n}}_{12,ij}$ and $\sigma^{2\hat{n}}_{12,ij}$, respectively. These relations will serve us below, since moments are more directly computed from the theory than cumulants.

Using eq. (3.9), the correlation function in redshift space is

$\displaystyle 1+\xi_{s}(\bm{s})$

$\displaystyle=\int d^{3}r\big{[}1+\xi(r)\big{]}\int\frac{d^{3}k}{(2\pi)^{3}}\exp\Big{[}i\bm{k}\cdot(\bm{s}-\bm{r}-\bm{v}_{12}^{\hat{\bm{n}}})-\frac{1}{2}\bm{k}^{T}\bm{\sigma}^{2\hat{\bm{n}}}_{12}\,\bm{k}+\cdots\Big{]}.$

(3.14)

If we stop at the second order cumulant $\bm{\sigma}^{2\hat{\bm{n}}}_{12}$, the $k$-integral can be formally performed analytically, giving

$$1+\xi_{s}(\bm{s})=\int\frac{d^{3}r}{(2\pi)^{3/2}|\mathbf{\sigma}^{2}_{12}|^{1/2}}\big{[}1+\xi(r)\big{]}\exp\left[-\frac{1}{2}(\bm{s}-\bm{r}-\bm{v}_{12}^{\hat{\bm{n}}})^{\mathbf{T}}[\bm{\sigma}^{2\hat{\bm{n}}}_{12}]^{-1}(\bm{s}-\bm{r}-\bm{v}_{12}^{\hat{\bm{n}}})\right],$$

(3.15)

which is the Gaussian Streaming Model correlation function.

Now, depending on how one computes the ingredients $\xi(r)$, $v_{12,i}(\bm{r})$ and $\sigma_{12,ij}^{2}(\bm{r})$, different methods can be adopted from here. For example: 1) Reference [40] computed $\xi(r)$ within the Zeldovich approximation, but the pairwise velocity ($v_{12}$) and pairwise velocity dispersion ($\sigma_{12}^{2}$) using Eulerian linear theory. 2) In [41] Convolution Lagrangian Perturbation Theory (CLPT) is used for the three ingredients, but instead of computing $\sigma_{12}^{2}$, the authors computed $\hat{\sigma}_{12}^{2}$. This latter reference also released a widely used code by the community.⁶⁶6Available at github.com/wll745881210/CLPT_GSRSD.

Here, we will use the method of [55, 56, 57], where all moments are computed using CLPT. Further, in our modeling we consider a Lagrangian biasing function $F$ that relates the galaxies and matter fluctuations through [58, 59]

$$1+\delta_{X}(\bm{q})=F(\delta,\nabla^{2}\delta)=\int\frac{d^{2}\mathbf{\Lambda}}{(2\pi)^{2}}\tilde{F}(\mathbf{\Lambda})e^{i\mathbf{D}\cdot\mathbf{\Lambda}}.$$

(3.16)

In the second equality $\tilde{F}(\mathbf{\Lambda})$ is the Fourier transform of $F(\bm{D})$, with arguments $\mathbf{D}=(\delta,\nabla^{2}\delta)$ and spectral parameters $\mathbf{\Lambda}=(\lambda,\eta)$, dual to $\mathbf{D}$. A key assumption that we follow here, is the number conservation of tracers,⁷⁷7Notice the number conservation assumption of tracers is not even true for halos. However, the biasing expansion obtained in this way is automatically renormalized and coincides with other more popular methods that introduce the biasing through the symmetries of the theory; see [60] for a review. from which one obtains

$$1+\delta_{X}(\bm{x})=\int\frac{d^{3}k}{(2\pi)^{3}}\int d^{3}qe^{i\bm{k}\cdot(\bm{x}-\bm{q})}\int\tilde{F}(\mathbf{\Lambda})e^{i\mathbf{D}\cdot\mathbf{\Lambda}-i\bm{k}\cdot\mathbf{\Psi}},$$

(3.17)

and evolves initially biased tracer densities using the map between Lagrangian and Eulerian coordinates given by eq. (3.1). Renormalized bias parameters are obtained as [58, 59]

$$b_{nm}=\int\frac{d\mathbf{\Lambda}}{(2\pi)^{2}}\tilde{F}(\mathbf{\Lambda})e^{-\frac{1}{2}\mathbf{\Lambda}^{\text{T}}\mathbf{\Sigma}\mathbf{\Lambda}}(i\lambda)^{n}(i\eta)^{m},$$

(3.18)

with covariance matrix components $\Sigma_{11}=\langle\delta_{cb}^{2}\rangle$, $\Sigma_{12}=\Sigma_{21}=\langle\delta\nabla^{2}\delta\rangle$ and $\Sigma_{22}=\langle(\nabla^{2}\delta)^{2}\rangle$. We notice $b_{n}=b_{n0}$ are local Lagrangian bias parameters, and $b_{\nabla^{2}\delta}=b_{01}$ is the curvature bias. In this work we consider only $b_{1}$ and $b_{2}$. However, tidal Lagrangian bias, $b_{s^{2}}$, can be easily introduced following [56].⁸⁸8Indeed, our code gsm-eft consider tidal bias, but we use the formulae presented in [61], which differ slightly from that in [56].

To obtain the cumulants we need to compute the moments. The procedure is exactly the same as with the correlation function (the zero order moment), but now we have to keep track of the velocity fields. That is, using

$$1+\delta_{X}(\bm{x}_{1})=\int d^{3}q_{1}\frac{d^{3}k_{1}}{(2\pi)^{3}}\frac{d\lambda_{1}}{2\pi}\tilde{F}(\lambda_{1})e^{i\lambda\delta_{1}+i\bm{k}_{1}\cdot(\bm{x}_{1}-q_{1}-\Psi_{1})},$$

we obtain

	$\displaystyle\Xi^{(n)}_{i_{1}\cdots i_{n}}(\bm{r})$	$\displaystyle=\hat{n}_{i_{i}}\cdots\hat{n}_{i_{n}}\int d^{3}q\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot(\bm{q}-\bm{r})}\int\frac{d\lambda_{1}}{2\pi}\frac{d\lambda_{2}}{2\pi}\tilde{F}(\lambda_{1})\tilde{F}(\lambda_{2})$
		$\displaystyle\quad\qquad\qquad\qquad\times\hat{n}_{j_{i}}\cdots\hat{n}_{j_{n}}\Big{\langle}\frac{\dot{\Delta}_{j_{1}}}{H}\cdots\frac{\dot{\Delta}_{j_{n}}}{H}e^{i[\lambda_{1}\delta_{1}+\lambda_{2}\delta_{2}+\bm{k}\cdot\Delta]}\Big{\rangle},$

where we defined $\Delta_{i}\equiv\Psi_{i}(\bm{q}_{2})-\Psi_{i}(\bm{q}_{1})$ and used $\Delta u_{i}=H^{-1}(\hat{n}_{j}\dot{\Delta}_{j})\hat{n}_{i}$.

The real space correlation function $\xi_{X}(r)$, which corresponds to the zeroth-order moment for tracer $X$, is obtained within CLPT, [62, 63, 12, 55, 64, 59],

		$\displaystyle 1+\xi_{X}(r)=\int\frac{d^{3}q}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}\Bigg{\{}1-\frac{1}{2}A_{ij}^{loop}G_{ij}+\frac{1}{6}\Gamma_{ijk}W_{ijk}$
		$\displaystyle\quad+b_{1}(-2U_{i}g_{i}-A^{10}_{ij}G_{ij})+b_{1}^{2}(\xi_{L}-U_{i}U_{j}G_{ij}-U_{i}^{11}g_{i})+b_{2}(\frac{1}{2}\xi_{L}^{2}-U_{i}^{20}g_{i}-U_{i}U_{j}G_{ij})$
		$\displaystyle\quad-2b_{1}b_{2}\xi_{L}U_{i}g_{i}+2(1+b_{1})b_{\nabla^{2}\delta}\nabla^{2}\xi_{L}+b^{2}_{\nabla^{2}\delta}\nabla^{4}\xi_{L}\Bigg{\}},$		(3.19)

where the matrix $A_{ij}(\bm{q})=\langle\Delta_{i}^{(1)}\Delta_{j}^{(1)}\rangle_{c}$, with $\Delta_{i}=\Psi_{i}(\bm{q}_{2})-\Psi_{i}(\bm{q}_{1})$, is the correlation of the difference of linear displacement fields for initial positions separated by a distance $\bm{q}=\bm{q}_{2}-\bm{q}_{1}$, which is further split in linear and loop pieces:

$$A_{ij}(\bm{q})=2\int\frac{d^{3}p}{(2\pi)^{3}}\big{(}1-e^{i\bm{p}\cdot\bm{q}}\big{)}\frac{p_{i}p_{j}}{p^{4}}P_{L}(p)=A^{L}_{ij}(\bm{q})+A^{loop}_{ij}(\bm{q}).$$

(3.20)

where $P_{L}$ is the linear matter power spectrum. We further use the linear (standard perturbation theory) correlation function

$$\xi_{L}(q)=\int\frac{d^{3}p}{(2\pi)^{3}}e^{i\bm{p}\cdot\bm{q}}P_{L}(p),$$

(3.21)

and the functions

$\displaystyle W_{ijk}=\langle\Delta_{i}\Delta_{i}\Delta_{k}\rangle_{c},\qquad A_{ij}^{mn}=\langle\delta^{m}(\bm{q})\delta^{n}(0)\Delta_{i}\Delta_{i}\rangle_{c},\qquad U^{mn}_{i}=\langle\delta^{m}(\bm{q})\delta^{n}(0)\Delta_{i}\rangle_{c}.$

(3.22)

The involved $r$ and $q$ dependent tensors are $g_{i}=(\mathbf{A}_{L}^{-1})_{ij}(r_{j}-q_{j})$, $G_{ij}=(\mathbf{A}_{L}^{-1})_{ij}-g_{i}g_{j}$, and $\Gamma_{ijk}=(\mathbf{A}_{L}^{-1})_{\{ij}g_{k\}}-g_{i}g_{j}g_{k}$.

The first and second moments of the generating function yield the pairwise velocity

		$\displaystyle v_{12,i}(\bm{r})=\frac{f}{1+\xi_{X}(r)}\int\frac{d^{3}q\,e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}\Bigg{\{}-g_{r}\dot{A}_{ri}-\frac{1}{2}G_{rs}\dot{W}_{rsi}$
		$\displaystyle\quad+b_{1}\left(2\dot{U}_{i}-2g_{r}\dot{A}^{10}_{ri}-2G_{rs}U_{r}\dot{A}_{si}\right)+b_{1}^{2}\left(\dot{U}^{11}_{i}-2g_{r}U_{r}\dot{U}_{i}-g_{r}\dot{A}_{ri}\xi_{L}\right)$
		$\displaystyle\quad+b_{2}\left(\dot{U}^{20}_{i}-2g_{r}U_{r}\dot{U}_{i}\right)+2b_{1}b_{2}\xi_{L}\dot{U}_{i}+2b_{\nabla^{2}\delta}\nabla_{i}\xi_{L}\,\Bigg{\}},$		(3.23)

and the pairwise velocity dispersion

		$\displaystyle\sigma^{2}_{12,ij}(\bm{r})=\frac{f^{2}}{1+\xi_{X}(r)}\int\frac{d^{3}q\,e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}\Bigg{\{}\ddot{A}_{ij}-g_{r}\ddot{W}_{rij}-G_{rs}\dot{A}_{ri}\dot{A}_{sj}$
		$\displaystyle\quad+2b_{1}\left(\ddot{A}^{10}_{ij}-g_{r}\dot{A}_{r\{i}\dot{U}_{j\}}-g_{r}U_{r}\ddot{A}_{ij}\right)+b_{1}^{2}\left(\xi_{L}\ddot{A}_{ij}+2\dot{U}_{i}\dot{U}_{j}\right)+2b_{2}\dot{U}_{i}\dot{U}_{j}\,\Bigg{\}},$		(3.24)

with $q$-coordinate dependent correlators

	$\displaystyle\dot{A}_{ij}^{mn}(\bm{q})=\frac{1}{fH}\langle\delta^{m}_{1}\delta^{n}_{2}\Delta_{i}\dot{\Delta}_{j}\rangle,\qquad\ddot{A}_{ij}^{mn}(\bm{q})=\frac{1}{f^{2}H^{2}}\langle\delta^{m}_{1}\delta^{n}_{2}\dot{\Delta}_{i}\dot{\Delta}_{j}\rangle,$
	$\displaystyle\dot{W}_{ijk}=\frac{1}{fH}\langle\Delta_{i}\Delta_{j}\dot{\Delta}_{k}\rangle,\qquad\ddot{W}_{ijk}=\frac{1}{f^{2}H^{2}}\langle\Delta_{i}\dot{\Delta}_{j}\dot{\Delta}_{k}\rangle,$
	$\displaystyle\dot{U}^{mn}(\bm{q})=\frac{1}{fH}\langle\delta^{m}_{1}\delta^{n}_{2}\dot{\Delta}_{i}\rangle.$		(3.25)

As in the case of the undotted $A_{ij}$ function, we have omitted to write the superscripts $m,n$ when these are zero; e.g, $\dot{A}_{ij}\equiv\dot{A}^{00}_{ij}$.

Now, consider the terms inside the curly brackets in eq. (3). Taking its large scales limit ($\bm{q}\rightarrow\infty$), we obtain

$\displaystyle\{\cdots\}\big{|}_{q\rightarrow\infty}$

$\displaystyle=\delta_{mn}\int_{0}^{\infty}\frac{dk}{3\pi^{2}}\left[P_{L}(k)+\text{loop terms}\right],$

(3.26)

which is a non-zero zero-lag correlator. However, since perturbation theory cannot model accurately null separations, one needs to add an EFT counterterm that has the same structure. This new contribution shifts the pairwise velocity dispersion as [56]

$$\hat{\sigma}^{2}_{12,mn}\rightarrow\hat{\sigma}^{2}_{12,mn}+\sigma^{2}_{\text{EFT}}\delta_{mn}\frac{1+\xi^{\text{ZA}}(r)}{1+\xi^{\text{CLPT}}_{X}(r)}.$$

(3.27)

As the separation distance $r$ increases, the ratio $(1+\xi^{\text{ZA}}(r))/(1+\xi_{X}(r))$ approaches unity, then the EFT counterterm adds as a constant shift to the pairwise velocity dispersion at large scales. That is, we can identify it with the phenomenological parameter $\sigma^{2}_{\text{FoG}}$ widely used in early literature to model Fingers of God (FoG). Comparisons for the modeling of the second moment when using the EFT parameter $\sigma^{2}_{\text{EFT}}$ and the constant shift $\sigma^{2}_{\text{FoG}}$ can be found in [61] (see for example Fig. 2 of that reference where a particular example exhibits a clear improvement of the EFT over the phenomenological constant shift). Finally, we notice our counterterm is related to that in [56] by $\sigma^{2}_{\text{EFT}}=\alpha_{\sigma}f^{2}$.

There are others EFT counterterms entering the CLPT correlation function and the pairwise velocity and velocity dispersion, but they are either degenerated with curvature bias (as is the case of $c_{1}^{\text{EFT}}$) or subdominant with respect to the contribution of eq. (3.27) (see the discussion in [56]). So, in this work we keep only $\sigma^{2}_{\text{EFT}}$.

Since this EFT parameter modifies the second cumulant of the pairwise velocity generation function, its effect on the redshift space monopole correlation function is small, while the quadrupole is quite sensitive to it, particularly at intermediate scales $r<40$ $h^{-1}\textrm{Mpc}$.

Now, let us comeback to eq. (3.14), that we formally integrated to obtain eq. (3.15). However, notice the matrix $\bm{\sigma}^{2\hat{\bm{n}}}_{12}$ is not invertible since $\bm{\sigma}^{2\hat{\bm{n}}}_{12}=\sigma^{2}_{12}\,\hat{\bm{n}}\otimes\hat{\bm{n}}$ and hence $\text{det}(\bm{\sigma}^{2\hat{\bm{n}}}_{12})=0$. Hence in the following we will approach this integration differently that will also serves us to rewrite the resulting equation in a more common form and also more directly related with the computational algorithms in a code.

We decompose the vectors $\bm{k}$, $\bm{s}$ and $\bm{r}$ in components parallel and perpendicular to the line of sight $\hat{\bm{n}}$:

$$\bm{k}=k_{\parallel}\hat{\bm{n}}+\bm{k}_{\perp},\quad\bm{k}=r_{\parallel}\hat{\bm{n}}+\bm{r}_{\perp},\qquad\bm{s}=s_{\parallel}\hat{\bm{n}}+\bm{s}_{\perp},$$

(3.28)

with $\bm{k}_{\perp}\cdot\hat{\bm{n}}=0$, and so on. We will use the following definitions

	$\displaystyle\mu$	$\displaystyle=\hat{\bm{r}}\cdot\hat{\bm{n}}=\frac{r_{\parallel}}{r},$		(3.29)
	$\displaystyle v_{12}(r)$	$\displaystyle=v_{12,i}\hat{r}_{i},$		(3.30)
	$\displaystyle\sigma^{2}_{12}(r,\mu)$	$\displaystyle=\mu^{2}\sigma^{2}_{12,\parallel}(r)+(1-\mu^{2})\sigma^{2}_{12,\perp}(r)$
		$\displaystyle=\mu^{2}(\hat{\sigma}^{2}_{12,\parallel}-v_{12}v_{12})+(1-\mu^{2})\hat{\sigma}^{2}_{12,\perp}(r),$		(3.31)

with $\hat{\sigma}^{2}_{12,\parallel}(r)=\hat{r}_{i}\hat{r}_{j}\hat{\sigma}^{2}_{12,ij}$ and $\hat{\sigma}^{2}_{12,\perp}(r)=\frac{1}{2}(\delta_{ij}-\hat{r}_{i}\hat{r}_{j})\hat{\sigma}^{2}_{12,ij}$, and

$$\bm{v}^{\hat{\bm{n}}}_{12}=\mu v_{12}(r)\hat{\bm{n}},\qquad\text{and}\qquad\bm{\sigma}^{2\hat{\bm{n}}}_{12}=\sigma^{2}_{12}(r,\mu)\hat{\bm{n}}\otimes\hat{\bm{n}},$$

(3.32)

Then, we can split the $\bm{k}$ integral eq. (3.14) in parallel and perpendicular to the line-of-sight integrations,

	$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot(\bm{s}-\bm{r}-\bm{v}^{\hat{\bm{n}}}_{12})-\frac{1}{2}\bm{k}^{T}\bm{\sigma}^{2\hat{\bm{n}}}_{12}\bm{k}}$	$\displaystyle=\int\frac{dk_{\parallel}}{2\pi}e^{ik_{\parallel}(s_{\parallel}-r_{\parallel}-\mu v_{12})-\frac{1}{2}k_{\parallel}^{2}\sigma^{2}_{12}}\int\frac{d^{2}k_{\perp}}{(2\pi)^{2}}e^{i\bm{k}_{\perp}\cdot(\bm{s}_{\perp}-\bm{r}_{\perp})}$
		$\displaystyle=\frac{e^{-\frac{1}{2\sigma^{2}_{12}}(s_{\parallel}-r_{\parallel}-\mu v_{12})^{2}}}{[2\pi\sigma^{2}_{12}]^{1/2}}\delta_{\text{D}}(\bm{s}_{\perp}-\bm{r}_{\perp}),$

obtaining a Dirac delta function from the integral of the perpendicular component $\bm{k}_{\perp}$ and a Gaussian kernel from the parallel $k_{\parallel}$ one.

Hence, the correlation function within the GSM becomes

$$1+\xi_{s}(s_{\parallel},s_{\perp})=\int_{-\infty}^{\infty}\frac{dr_{\parallel}}{\big{[}2\pi\sigma^{2}_{12}(r,\mu)\big{]}^{1/2}}\big{[}1+\xi(r)\big{]}\exp\left[-\frac{\left(s_{\parallel}-r_{\parallel}-\mu v_{12}(r)\right)^{2}}{2\sigma^{2}_{12}(r,\mu)}\right],$$

(3.33)

with $r^{2}=r_{\parallel}^{2}+s_{\perp}^{2}$. This is a wide popular expression, but remind that here $\sigma^{2}_{12}$ is the second cumulant of the density weighted velocity generating function, instead of its second moment. Also, it suffers correction from EFT counterterms.

The streaming models [65, 38, 39] describe how the fractional excess of pairs in redshift space $1+\xi_{s}(\bm{s})$ is modified with respect to their real-space counterpart $1+\xi(r)$:

$$1+\xi_{s}(s_{\parallel},s_{\perp})=\int_{-\infty}^{\infty}dr_{\parallel}\big{[}1+\xi(r)\big{]}\mathcal{P}(s_{\parallel}-r_{\parallel}|\bm{r}).$$

(3.34)

Here $r^{2}=r^{2}_{\parallel}+r^{2}_{\perp}$ and $s_{\perp}=r_{\perp}$. The above expression is exact; see eqs.(1)-(12) of [39]. This means that a knowledge of the form of the pairwise velocity distribution function $\mathcal{P}(v_{\parallel}|\bm{r})=\mathcal{P}(s_{\parallel}-r_{\parallel}|\bm{r})$ at any separation $\bm{r}$, yields a full mapping of real- to redshift-space correlations. In the GSM approximation, the distribution function becomes Gaussian centered at $\mu v_{12}$ and with width equal to $\sigma_{12}$.

The main drawback of this approach is that the $\mathcal{P}(v_{\parallel}|\bm{r})$ is, of course, not a Gaussian [39, 40, 66]. In [66], the authors extract the ingredients of the pairwise velocity disribution moments directly from simulations and use them to obtain the correlation function multipoles using the GSM and the Edgeworth streaming model of [55], finding good agreement with the redshift space correlation function extracted from the same simulations, but only above scales of around $s=20\,h^{-1}\text{Mpc}$. Our findings also indicate that our modeling and pipeline fits well the simulations above this same scale.

Our full shape analysis requires an exploration of a relatively large parameter space. Each model requires approximately 1.5 seconds in our computer to run. Given the large number of evaluations required to explore the parameter space (of the order of $10^{5}$), and the large amount of MCMC chains we are interested in running, there is an incentive to optimize the evaluation process of our model.

There are various methodologies available for accelerating the estimation of these statistics. The choice between using one or another depends on several factors, one of them being the number of models that can be constructed to use as a training set. In our case, the Gaussian streaming model presented in Section 3 is relatively cost-efficient. To run our model, we first construct a template of the power spectrum using the publicly available Code for Anisotropies in the Microwave Background (CAMB) [67, 68], which completes in approximately one second. We then utilize this template as input for our gsm-eft code, which requires an additional half-second to compute the correlation function multipoles. Considering this, we can efficiently generate training data sets of several tens of thousands of points within a reasonable computational time. Recently, neural networks, have proved to be a suitable framework to accelerate the estimation of clustering statistics for training sets of this size [e.g. 69, 37]. Moreover, neural networks are particularly efficient in data generalization, affording an almost constant reliability over the full parameter space (i.e. the model error does not strongly depends on the distance with the nearest point used on the training set).

In what follows, we present our emulating methodology. Our approach is derived from the methodology proposed in [37], but we have made modifications to adapt it for configuration space. The following subsection provides a detailed explanation of our methodology and highlights the specific changes we have implemented to transition into configuration space. Once our neural networks are trained, we reduce the evaluation time needed for a single point in parameter space to around 0.015 seconds, which improves in two orders of magnitude the Likelihood evaluation time.

We note that a distinct neural network emulator is necessary for each multipole of the correlation function at a specific redshift.¹¹¹¹11One can decide to train a global neural network including the two multipoles. However, it corresponds to expanding the exit layer by a factor of three without a win of information between them. Throughout this study, we utilize the first two non-zero multipoles of the correlation function. Consequently, each analysis presented in this work entails constructing two neural network emulators. The construction process for each emulator takes approximately 30 minutes when performed on our personal laptops. It is worth mentioning that it might be feasible to reduce the building time of a single neural network by adjusting certain parameters, as discussed below. Nevertheless, since the number of neural networks we use is small, we find the current building time to be manageable.

Other methodologies that operate in Fourier space [e.g. 70, 37], aim to predict values for all wave numbers of interest, which typically comprise hundreds of points. If a brute force approach were employed, where one asks the neural network to directly predict the power spectrum, it would need to make hundreds of predictions, which would increase the time required to build the network. Therefore, methodologies utilizing Fourier space often employ techniques such as principal component analysis to address this issue.¹²¹²12In principal component analysis, the input power spectra matrix is divided into eigenvectors, which are dependent on the wave number, and their corresponding eigenvalues, which only rely on the cosmology. This enables an approximation of the power spectra by considering a linear combination of the most significant eigenvalues and discarding the rest, reducing the number of predictions necessary.

Here, we model the correlation function from 20 $h^{-1}\,\textrm{Mpc}$ to 130 $h^{-1}\,\textrm{Mpc}$ in 22 bins with a 5 $h^{-1}\,\textrm{Mpc}$ width between each bin in redshift space distance $r$, therefore, each emulator needs to predict only 22 numbers.

We have found no necessity to incorporate principal component analysis as part of our methodology as Fourier space works utilize a number of principal components similar to the number of bins we employ.

To train each neural network, we generate 60,000 models distributed across the parameter space. Out of these, 50,000 models are utilized for training, while 5,000 models constitute the validation set. The validation set is employed during the training process to test the data and determine when to decrease the learning rate of the network, as explained below. The remaining 5,000 models form our test set and are reserved for evaluating the accuracy of the methodology on unseen data so that we perform a fair assessment of the trained models.

We use Korobov sequences [71] to select the points in the parameter space used for building and testing our neural networks. Korobov sequences are a robust approach for generating extensive and uniform samples in large dimensional spaces. We run three distinct sequences to create the training, test, and validation sets at each distinct redshift that we model. We also make sure that the three sequences are independent and that there are no overlapping points between them. Finally, we employ our EFT-GSM model pressented in section 3 to calculate the multipoles of the power spectra for all 60,000 data points. These EFT-GSM multipoles are used to train our neural networks, that aim to accurately replicate the GSM predictions for a new point in parameter space.

In order to make the training of the neural network more efficient, it is convenient to keep the values of the output layer neurons in a similar range of values. Here, we use a hyperbolic sinus transformation on each of the training set multipoles for this purpose.

We run our neural networks by adapting the public code from [37],¹³¹³13Which is available at https://github.com/sfschen/EmulateLSS/tree/main to reflect the changes expressed above. We use the Multi-Layer Perceptron architecture suggested by them, with four hidden layers of 128 neurons each. As we discuss below, the accuracy that we obtain in our predictions is well within the precision we need and this is achieved in a manageable time. Therefore, we decided that no further optimization of the architecture to fit our particular data was necessary. When training our networks, we reduce the learning rate of the algorithm from $10^{-2}$ to $10^{-6}$ in steps of one order of magnitude and double the training batch size at every step. As suggested by several works [e.g. 72, 37], we use the following activation function,

$$a(X)=\left[\gamma+(1+e^{-\beta\odot X})^{-1}(1-\gamma)\right]\odot X,$$

(4.1)

which we found outperforms other more common activation functions like Rectified linear units. Here, $\gamma$ and $\beta$ are new free parameters of a given hidden layer within the neural network that are fitted during the training process of the network.

The algorithm decreases the learning rate when a predetermined number of training epochs have passed without any significant improvement on the accuracy of the model, this number is commonly referred to as the patience of the algorithm.¹⁴¹⁴14We track the mean square error (MSE) of the validation set, the algorithm records the best value found so far. When a number of epochs equal to our patience value have elapsed whiteout a better MSE being found, the algorithm switches the learning rate. Note that a larger patience allows the model more time to exit local minima, and address slow convergence issues. The patience we use determines the time required to train our neural networks, longer patience usually leads to more accurate models (provided the model is not overfitted). The results presented in this work correspond to waiting 1000 epochs before reducing the learning rate, which, as stated above corresponds to approximately 30 minutes of training time, also, we have monitored our validation set to ensure that there are no signs of overfitting at this point.¹⁵¹⁵15 An overfitted model would start to worsen the MSE of the validation set after a given training epoch. If our goal were to reduce the training time of the algorithm, we could decrease the patience value. However, this would reduce the accuracy of our models, we note that a patience of around 100 epochs reduces the training time to around 5 minutes on our personal laptops while still maintaining sub-percent accuracy in most multipole models.

In this section, we define the methodology employed to extract cosmological information from galaxy clustering. First, we describe the clustering measurements we utilize. Next, we provide an overview of our full shape methodology. We also discuss the likelihood, priors, covariance, and MCMC samplers used throughout our analysis.