Cosmological constraints from the full-shape galaxy power spectrum in SDSS-III BOSS using the BACCO hybrid Lagrangian bias emulator Open Access

ABSTRACT

We present a novel analysis of the redshift-space power spectrum of galaxies in the Sloan Digital Sky Survey III (SDSS-III) Baryon Oscillation Spectroscopic Survey (BOSS). Our methodology improves upon previous analyses by using a theoretical model based on cosmological simulations coupled with a perturbative description of the galaxy–matter connection and a phenomenological prescription of fingers of God. This enables a very robust analysis down to mildly non-linear scales, |$k\simeq 0.4 \, h\, {\rm Mpc}^{-1}$|⁠. We carried out a number of tests on mock data, different subsets of BOSS, and using model variations, all of which support the robustness of our analysis. Our results provide constraints on |$\sigma _8$|⁠, |$\Omega _m$|⁠, h, and |$S_8 \equiv \sigma _8 \sqrt{ \Omega _{\rm m} /0.3}$|⁠. Specifically, we measure |$\Omega _m=0.301\pm 0.011$|⁠, |$\sigma _8=0.745^{+0.028}_{-0.035}$|⁠, |$h=0.705\pm 0.015$|⁠, and |$S_8 = 0.747^{+0.032}_{-0.039}$| when all the nuisance parameters of our model are left free. By adopting relationships among bias parameters measured in galaxy formation simulations, the value of |$S_8$| remains consistent whereas uncertainties are reduced by |$\sim 20~{{\ \rm per\ cent}}$|⁠. Our cosmological constraints are some of the strongest obtained with the BOSS power spectrum alone: they exhibit a |$2.5{\!-\!}3.5\sigma$| tension with the results of the Planck satellite, agreeing with the lower values of |$S_8$| derived from gravitational lensing. However, the cosmological model preferred by Planck is still a good fit to the BOSS data, assuming small departures from physical bias priors and, therefore, cannot be excluded at high significance. We conclude that, at the present, the BOSS data alone does not show strong evidence for a tension between the predictions of Lambda cold dark matter (⁠|$\Lambda$|CDM) for the high- and low-redshift Universe.

1 INTRODUCTION

The spatial distribution of galaxies encodes information about fundamental aspects of our Universe. For instance, the clustering amplitude as a function of scale depends on the abundance of different elements – dark matter, baryons, and photons – as well as the properties of the primeval Universe. Additionally, the clustering anisotropy induced by the peculiar velocities of galaxies, known as redshift-space distortions (RSDs), depends on the growth rate of cosmic structures, which contains information about cosmological parameters and the law of gravity.

This rich source of information has motivated large spectroscopic surveys, which map the three-dimensional distribution of up to hundreds of millions of galaxies. It is not an exaggeration to say that a revolution in cosmology is on the horizon. Recent images released from ESA’s EUCLID satellite (Euclid Collaboration 2024) show the ability of the mission to deliver unprecedented maps of the cosmos. Additionally, after only 1 yr of observations, the Dark Energy Spectroscopic Instrument (DESI; Levi et al. 2013) already built the largest ever galaxy redshift catalogue. Over the next decade, these surveys will be complemented by a myriad of other experiments such as 4MOST (de Jong et al. 2019), SPHEREX (Crill et al. 2020), and J-PAS (Bonoli et al. 2021).

An essential aspect of a successful exploitation of clustering data is accurate theoretical models. On the one hand, the spatial distribution of galaxies can be described with analytical approaches. These models are based on perturbative solutions of the relevant evolution equation for matter and for the connection between galaxies and the matter field (see the reviews by Bernardeau et al. 2002 and Desjacques, Jeong & Schmidt 2018). Although developments, such as the effective field theory (EFT) of Large Scale Structure (LSS; e.g. Baumann et al. 2012; Carrasco et al. 2014; Schmidt 2021; Cabass et al. 2023), have considerably improved the accuracy and range of validity of such models, these approaches are only valid on large scales. On the other hand, in recent years it is becoming feasible to make direct use of cosmological simulations in performing analyses of redshift-space clustering data. These approaches have mostly been based on assuming Halo Occupation Distribution models, where the abundance of galaxies is primarily given by the mass of the host halo (see e.g. Sullivan, Seljak & Singh 2021; Hahn et al. 2022; Cuesta-Lazaro et al. 2023; Yuan, Abel & Wechsler 2024). Although these models can predict the clustering down to much smaller scales than, e.g. EFT and therefore are potentially much more constraining, they make rigid assumptions about galaxy formation physics, which might not apply to the real universe. This could bias cosmological inferences, especially for samples with complex selection criteria.

An alternative family of theoretical models is given by hybrid methods (Modi, Chen & White 2020), which combine the accuracy of numerical simulations with the general formulation of pertubative models. Specifically, perturbative approaches are used to describe in full generality the large-scale connection between galaxies and the linear matter field in Lagrangian coordinates. This Lagrangian ‘biased’ field is then mapped to Eulerian space, and subsequently to redshift space (Pellejero Ibañez et al. 2022), with N-body simulations. Therefore, we can take advantage of the accuracy of numerical simulations to increase the performance of perturbative models. Recently, this model has been used in studies of the |$k-$|NN statistics in Banerjee, Kokron & Abel (2022), and 21-cm field in Baradaran et al. (2024).

In this paper, we will employ the ‘BACCO hybrid Lagragian bias’ emulator that implements a second-order Lagrangian bias expansion on top of the BACCO simulation suite (Angulo et al. 2021). We have developed and validated our emulator in several previous papers. In Zennaro et al. (2021) (see also Kokron et al. 2021 for an independent confirmation for clustering and lensing), we showed that our approach describes the clustering of mock galaxies down to quasi-nonlinear scales. Then, in Zennaro et al. (2022) we validated it against thousands of physical galaxy mocks built with various flavours of subhalo abundance matching. In parallel, in Pellejero Ibañez et al. (2022) we extended and validated the hybrid-bias model for redshift-space clustering, and in Pellejero Ibañez et al. (2024) we explored its combination with field-level emulators. Finally, in Pellejero Ibañez et al. (2023) we built the redshift-space model for thousands of different sets of cosmological parameters that included massive neutrinos and dynamical dark energy. With such data, we trained a neural-network emulator that predicts the monopole, quadrupole, and hexadecapole of the power spectrum for biased tracers as a function of cosmology. Our emulators have been further tested and compared with other approaches in preparation for the EUCLID and Rubin LSST surveys (Euclid Collaboration 2023; Nicola et al. 2024), and in the ‘Beyond-2point’ challenge (Beyond-2pt Collaboration 2024). An independent application to the DES Y1 data was done by Hadzhiyska et al. (2021).

The goal of this paper is to present the first application of our redshift-space emulator to observational data. For this, we will use public measurements of the monopole, quadrupole, and hexadecapole of the redshift-space power spectrum of galaxies measured in Sloan Digital Sky Survey III (SDSS-III) Baryon Oscillation Spectroscopic Survey (BOSS; Alam et al. 2017). This is the largest publicly available redshift galaxy survey, with a volume of nearly 6 |$[h^{-1}{\rm Gpc}]^3$| and over one million galaxies.

We will show that our results provide cosmological constraints in broad agreement with previous analyses using the same BOSS datavector and either perturbation theory or simulation-based models. On the other hand, our analysis delivers 30  per cent tighter constraints on |$\sigma _8$| with respect to pertubative approaches, which is a consequence of the unique features of our modelling.

Although our constraints on |$\Omega _m$| and |$H_0$| are consistent with the recent analysis of the Baryon Acoustic Oscillations (BAO) feature in the first year of data in DESI (DESI Collaboration 2024), our marginalized constraints on |$S_8$| differ by |$2.5-3.5\sigma$| with respect to those from Planck. Likewise, our |$S_8$| constraint is in tension with the results from Cosmic Microwave Background (CMB) lensing and cluster counts performed by the Atacama Cosmology Telescope (ACT, Madhavacheril et al. 2024) and the extended ROentgen Survey with an Imaging Telescope Array (eROSITA, Ghirardini et al. 2024) collaborations, respectively. However, our measurements are consistent with other low-redshift cosmological probes, such as joint analyses of lensing data and cluster abundances from the South Pole Telescope survey (SPT, Bocquet et al. 2024), other full-shape analysis (Philcox & Ivanov 2022; Chen, Vlah & White 2022a; Ivanov et al. 2023), and joint analyses of lensing data and galaxy clustering (Chen et al. 2022b, 2024a). Despite this tension, we note that the cosmological parameter set preferred by Planck represents a very good description of the BOSS data, which suggests that the current BOSS data do not show significant enough evidence for a fundamental tension with the high-z Universe within the standard Lambda cold dark matter (⁠|$\Lambda$|CDM) model.

Our paper is organized as follows. We start in Section 2 by briefly describing the characteristics of the BOSS measurement, and, in Section 3, recapping our model. In Section 4, we then validate our inference pipeline against a suite of N-body mock catalogues resembling the BOSS data. Subsequently, in Section 5 we apply our pipeline to the redshift-space power spectrum of BOSS galaxies. We will first consider wide priors on the relevant cosmological, bias, and nuisance parameters. In a second stage, we will perform an exploratory analysis where we fix the relationship between the linear and high-order bias parameters according to the expectations from galaxy simulations. In Section 6, we present several test of the robustness our results. Finally, we discuss the implications of our findings in Section 7, and in Section 8 we conclude.

2 GALAXY CLUSTERING DATA

In this paper, we will employ the redshift-space power spectrum of galaxies in SDSS-III BOSS. Below we provide a summary of the survey and the estimation of the power spectrum and the corresponding covariance matrix.

2.1 Baryon Oscillation Spectroscopic Survey

The data set corresponds to the |$12{\rm th}$| data release (DR12) of the BOSS (Dawson et al. 2013; Reid et al. 2016), a part of SDSS-III (Eisenstein et al. 2011). BOSS encompasses 1198 006 galaxies observed across two distinct regions of the sky: the Northern and Southern galactic caps (referred to as NGC and SGC, respectively). There are two galaxy selection criteria in BOSS: CMASS and LOWZ, which roughly correspond to high-redshift galaxies above a stellar mass threshold and low-redshift galaxies above a magnitude limit.

For the purposes of this analyses, we consider a division of each galactic cap only in terms of redshift: |$z_1 \in [0.2, 0.5]$| and |$z_3 \in [0.5, 0.75]$|⁠, with a median redshift of 0.38 and 0.61, respectively. Although each of these redshift ranges will contain a mixture of CMASS and LOWZ galaxies, we expect this not to affect our modelling given the generality of the bias expansion. The raw volume of each of these four subsamples is 1.7 and 3.3 Gpc|$^3$| for SGC, and 4.7 Gpc|$^3$| and 9 Gpc|$^3$| for the NGC (Alam et al. 2017), to give a total of 18.7 Gpc|$^3$| (approximately 6 |$[h^{-1}{\rm Gpc}]^3$|⁠). Therefore, we expect the high-redshift NGC to be the sample with the largest statistical power. Conversely, the low-redshift SGC is expected to display the weakest constraining power.

2.2 Power-spectrum estimation

We employ the power spectrum of galaxies in each of the four samples, as measured by Philcox & Ivanov (2022). The statistic is computed using an ‘unwindowed’ estimator, which aims to deliver the power spectrum without the effect of the survey mask, easing the comparison with theoretical models. In practice, the estimators are obtained by defining the large-scale likelihood for the galaxy survey and subsequently optimizing analytically for the desired statistics, analogous to methodologies employed in early analyses of CMB and LSS.

The specific formulation of the estimator used by Philcox & Ivanov (2022) was proposed in Philcox (2021a, b) and Philcox & Flöss (2024). These authors approximate the pixel covariance using a diagonal (FKP-like) structure with |$P_{\rm fkp} \approx 10^4\, [h^{-1}{\rm Mpc}]^3$|⁠, while still fully considering the survey’s geometry. In addition, the measurements are corrected by the impact of observational systematic errors such as galaxy-star separation, redshift failures, spectroscopic fibre collisions, and observing conditions. The public measurements¹ provide the monopole, quadrople, and hexadecapole of the power spectrum with a resolution of |$\Delta k = 0.005 \, h\, {\rm Mpc}^{-1}$| up to |$k=0.41 \, h\, {\rm Mpc}^{-1}$| which is given by the Nyquist frequency of the Fourier grid employed |$k_{\rm Nyq} = 0.45 \, h\, {\rm Mpc}^{-1}$|⁠. The Poisson contribution coming from the finite galaxy sample is removed from the observed statistics, which will impact the results of the noise contribution, as we will see later.

Note that the smallest scale provided by the public measurements will determine the wavenumbers we employ in our analysis, despite our theoretical model being, in principle, able to use even smaller scales (⁠|$k_{\rm max} \sim 0.6\, \, h\, {\rm Mpc}^{-1}$|⁠), as shown in Pellejero Ibañez et al. (2023).

2.3 Covariance matrix

The covariance matrix of the power-spectrum measurements was estimated using a set of 2048 multidark-patchy mock catalogues (Kitaura et al. 2016; Rodríguez-Torres et al. 2016), publicly provided by the BOSS collaboration.² In particular, we use the window-free estimation by Philcox & Ivanov (2022).

Each Patchy catalogue is built by creating a non-linear density and velocity field with a fast but approximate gravity solver, which is then calibrated to reproduce the results from a high-resolution N-body simulation. The Patchy catalogues have a flexible model for the galaxy–matter connection, whose free parameters were set to reproduce the observed redshift-space power spectrum in BOSS. Finally, mock galaxies are selected to display the same observational set-up as BOSS in terms of its window mask and weights due fibre collisions.

We note that the Patchy mocks have only been extensively validated for |$k\lt 0.3 \, h\, {\rm Mpc}^{-1}$|⁠, yet we will employ them until |$k \simeq 0.4 h^{-1}{\rm Mpc}$|⁠. We justify this by noting that, for the number density and volume of BOSS, the diagonal elements of the covariance matrix are dominated by the Gaussian terms, whereas off-diagonal elements are dominated by mode coupling caused by the survey mask (Colavincenzo et al. 2018; Lippich et al. 2018; Blot et al. 2019). Therefore, since all these elements have been carefully incorporated into the MultiDark-Patchy mocks, we expect our covariance matrix to provide a reasonably accurate description of the uncertainties over the full range of scales we will employ. In the future, we will study the dependency of our results with alternative methods to estimate covariance matrices (Balaguera-Antolínez et al. 2019a, b; Pellejero-Ibañez et al. 2020a; Sinigaglia et al. 2021; Kitaura et al. 2022; Ereza et al. 2023).

Although Patchy represents an adequate method to estimate covariance matrices, its inherent approximations might compromise the accuracy of the mock multipoles, specially at |$k \gt 0.2 \, h\, {\rm Mpc}^{-1}$|⁠. For this reason, we will validate our inference pipeline using a different suite of mock catalogues (see Section 4).

3 THEORETICAL MODELLING

In this section, we provide details about our cosmology inference pipeline. In Section 3.1, we discuss our modelling for the redshift-space clustering of biased tracers, including the hybrid bias expansion and how we emulate our predictions. In Section 3.2, we provide details on our Bayesian inference, including the definition of the parameter space, the form of the likelihood, and the sampler.

3.1 The galaxy power spectrum in redshift space

3.1.1 Hybrid Lagrangian bias expansion

The modelling of clustering statistics within the ‘hybrid’ approaches, as outlined in Modi et al. (2020) and extended in Pellejero Ibañez et al. (2022), involves two key components: (i) a mapping from Lagrangian (⁠|$\pmb {q}$|⁠) to Eulerian (⁠|$\pmb {x}$|⁠) space given by the displacement field, |$\pmb {\psi }(\pmb {q})$|⁠, measured in N-body simulations, and (ii) a functional relationship in Lagrangian space between the linear overdensity field, |$\delta _{\rm L}(\pmb {q})$|⁠, and the galaxy density field, |$\delta _{\rm g} (\pmb {q})$|⁠.

The relationship between matter and galaxies, |$F[ \delta _{\rm L}(\pmb {q}) ]$|⁠, is typically referred to as the bias model. In our approach, we adopt a second-order perturbative expansion for |$F(\pmb {q})$|⁠, as described in Matsubara (2008) and Desjacques et al. (2018), which weighs the significance of various Lagrangian fields in depicting the density of tracers at a given |$\pmb {q}$|⁠:

Here, brackets |$\langle \cdot \rangle$| denote volume averages, and |$s^2 \equiv s_{ij}s^{ij}$| where |$s_{ij}$| is the traceless tidal tensor. Thus, |$s^2 = \partial _i\partial _j\phi (\pmb {q})-1/3\delta ^{\rm {K}}_{ij}\delta _{\rm {L}}(\pmb {q})^2$|⁠, with |$\phi (\pmb {q})$| representing the linear gravitational potential.

Therefore, the overdensity field of galaxies in Eulerian space, |$\delta _{\rm g} (\pmb {x})$|⁠, reads

where |$\delta _{\rm D}$| is a Dirac’s delta. We highlight that our implementation of the hybrid bias expansion has undergone extensive testing. For instance, in Zennaro et al. (2022), we built thousands of mock catalogues based on the SubHalo Abundance Matching extended model (SHAMe Contreras, Angulo & Zennaro 2021a, b) with different parameter sets, selection criteria, number densities, redshifts, and cosmologies. We then demonstrated that our approach was able to accurately fit the galaxy power spectra down to |$k_{\rm max} \simeq 0.7 \, h\, {\rm Mpc}^{-1}$|⁠. Given the success of the model in real space, recently we extended the hybrid approach to model the intrinsic alignment of galaxy shapes (Maion et al. 2023).

To compute the galaxy overdensity in redshift space, |$\delta _{\rm g} ^z$|⁠, we include an additional term in the displacement field to account for the effect of peculiar velocities (Pellejero Ibañez et al. 2022):

where |$\hat{ \pmb {x} }_z({ \pmb {q} })$| represents the unit vector along the line of sight, a denotes the scale factor, and H is the Hubble parameter. The velocity field |$\pmb {v}$| is built directly from the peculiar velocities measured for haloes and matter in N-body simulations. Thus, our model inherits in full the non-linearity of velocities and of the mapping between real-, |$\pmb {x}$|⁠, and redshift-space coordinates, |$\pmb {s}$|⁠.

We note that we cannot simply use |$v(\pmb {x})$| as the velocity of matter in simulations because we expect galaxies to sample a coarsed-grained version of the velocity field. For instance, central galaxies will be mostly at rest relative to their host halo, whereas satellite galaxies could display larger or smaller velocity dispersion than typical particles in a halo, owing to the details of galaxy formation physics (Orsi & Angulo 2018; Alam et al. 2021). Therefore, we define |$\pmb {v}(\pmb {x}) = \pmb {v}_{\rm matter}(\pmb {x})$| if |$\pmb {x}$| is outside a halo and |$\pmb {v}(\pmb {x}) = \pmb {v}_{\rm {halo}}(\pmb {x})$| if |$\pmb {x}$| is found within a halo. To model the role of intra-halo satellite velocities, we include two additional free parameters that describe the fraction of satellite galaxies, |$f_{\rm sat}$|⁠, and their typical velocity dispersion |$\lambda _{\rm sat}$|⁠:

where |$\boldsymbol{\ast }_z$| denotes a convolution along the line-of-sight direction. In Pellejero Ibañez et al. (2022), we extensively tested this model for different SHAMe catalogues in redshift space (including one that mimics the CMASS sample), finding it describes the clustering of galaxies down to |$k_{\rm max} \sim 0.6 \, h\, {\rm Mpc}^{-1}$| at the precision of the BOSS survey. We also tested this model in a blind manner within the ‘Beyond 2-pt Challenge’ (Beyond-2pt Collaboration 2024), where our approach delivered unbiased and among the strongest constraints on |$\sigma _8$|⁠.

Although the model provides predictions at the field level (Pellejero Ibañez et al. 2024), in this paper we are only interested in the redshift-space power spectrum which is typically characterized by a multipole expansion:

where |$k\equiv |\pmb {k}|$|⁠, |$\mu \equiv \pmb {k} \cdot \pmb {k_z}$|⁠, and |$\mathcal {P}_\ell (\mu)$| are Legendre polynomials.

Finally, we include two stochastic terms in the power spectrum monopole, |$\ell =0$|⁠, to account for discreteness noise and small-scale physics not included in our bias model: |$\epsilon (k) = 1 /\bar{n} \, (\epsilon _1+\epsilon _2k^2)$|⁠, where |$\bar{n}$| is the mean number density of galaxies in the sample. Note that in Pellejero Ibañez et al. (2022) we tested the need to include |$\mu$|-dependent stochastic terms, finding them not required for the range of scales and level of accuracy of BOSS. This does not imply the term is absent, rather, it indicates that the velocity from haloes in the simulation, combined with the FoG free parameters, accounts for this dependency.

3.1.2 Emulation

For our analysis, we will employ the redshift-space hybrid bias emulator presented in Pellejero Ibañez et al. (2023). Below, we will recap its main features.

For our model to make predictions as a function of cosmological parameters, we require displacement and velocity fields measured in N-body simulations for multiple cosmologies. Obtaining such predictions is a computationally challenging task, thus, we resort to the BACCO suite of high-resolution simulations (initially introduced in Angulo et al. 2021), which can be transformed to span thousands of different cosmologies by employ a cosmology-rescaling algorithm (Angulo & White 2010; Zennaro et al. 2019; Contreras et al. 2020).

We highlight that, in addition to the large volume of BACCO simulations (⁠|$V\simeq 8 {\rm Gpc}^3$| with |$4320^3$| particles of |$\sim 10^{9} h^{-1}{\rm M_{ \odot }}$| mass), these simulations have ‘fixed-and-paired’ initial conditions, which significantly reduced the cosmic variance errors in the simulated displacement and velocity fields (Angulo & Pontzen 2016; Maion, Angulo & Zennaro 2022).

We rescaled the BACCO suite to 4000 combinations of cosmological parameters and redshifts, over a hyperspace in the cosmological parameters |$\lbrace \Omega _{\rm m}, \sigma _8 , \Omega _b, n_s, h, M_{\nu } \, [{\rm eV}], w_{0}, w_{a}\rbrace$| that roughly spans 10 times the uncertainty delivered by the analysis of Planck data (Planck Collaboration VI 2,020). Although our emulator provides predictions as a function of neutrino mass, |$M_{\nu }$|⁠, and dynamical dark energy, |$w(a) \equiv w_0 + (a-1)\, w_a$|⁠, we will perform our current analysis within the minimal |$\Lambda$|CDM model: |$M_{\nu } = 0$|eV, |$w_0=-1$|⁠, and |$w_a=0$|⁠.³

At each parameter set, we build the hybrid model predictions by advecting each of the 5 terms in equation (1) to redshift space, and compute all the 15 cross-power spectra that enter in equation (5). To reduce statistical noise as much as possible, we repeated and average our calculations over the three line-of-sight directions, |$\lbrace \hat{x}, \hat{y}, \hat{z}\rbrace$|⁠.

These measurements provided a library of power spectrum multipoles at different cosmologies that allows for emulation via a neural network. Pellejero Ibañez et al. (2023) demonstrated that the typical accuracy for the monopole, quadrupole, and hexadecapole is 0.5 per cent, 1 per cent, and 10 per cent, respectively. These are subdominant compared to the statistical uncertainties in BOSS, where the smallest uncertainties for the NGC subsample are |$3{\!-\!}4~{{\ \rm per\ cent}}$| for the monopole, 20 per cent for the quadrupole, and 40 per cent for the hexadecapole.

3.1.3 The Alcock–Paczynski effect

The clustering of galaxies would appear anisotropic if the cosmology assumed to transform redshifts and angles into comoving separations is incorrect.

To account for this effect, known as Alcock–Paczynski distortion (Alcock & Paczynski 1979), we first construct a two-dimensional power spectrum, |$P(k,\mu)$|⁠, from the three multipoles that we emulate. Then, we transform our measurements to how they would have appeared if we had assumed the same cosmology as the BOSS collaboration (Ballinger, Peacock & Heavens 1996):

where H and |$D_A$| are the Hubble parameter and the angular diameter distance at the median redshift of each sample. The ‘fid’ subscript refers to quantities computed adopting the fiducial cosmology in the BOSS survey.

At each likelihood evaluation, we compute the transformation at the corresponding cosmological parameters. We then compute |$P_{\rm fid}$| and the respective power-spectrum multipoles using equation (5).

3.2 Parameter inference

3.2.1 Likelihood and sampling

Due to the central limit theorem, we will assume that the likelihood of observing a set of power spectrum multipoles, |$\vec{P}_{\ell }=\lbrace P_0,P_2,P_4\rbrace$|⁠, has a multivariate Gaussian form:

where |$\vec{P}_\ell ^{\rm bacco}$| are the multipoles evaluated with the model described in Section 3.1.1, and |$C_{\rm {BOSS}}$| is the covariance matrix estimated using Patchy mocks (Kitaura et al. 2016), as detailed in Section 2.3. When combining samples, we assume them to be independent from each other, i.e. the full covariance matrix is block diagonal.

We compute the likelihood and credibility intervals on our parameters using a nested sampler algorithm, as implemented in the publicly available code multinest⁴ (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2019). Our analysis uses an evidence tolerance of 0.1 and a number of live points of 250. We tested with more stringent set-ups (live points of 400) for our fiducial result and found no differences, as opposed to Lemos et al. (2023). This is likely due to the reduced parameter space in comparison to their work.

After evaluating the 15 emulators, including the FoG and noise effects and performing the AP transformations, one likelihood evaluation takes approximately 1 s. This is a reasonable time for running chains with the nested sampler, provided we restrict ourselves to parameter spaces smaller than about 30 parameters. For larger parameter spaces, we would need to speed up the likelihood computation, for example, by including the AP parameters as part of the emulated parameters. We will explore these options in future work.

3.2.2 Parameter space

Our theoretical model, |$P_\ell ^{\rm bacco}$|⁠, contains 13 free parameters: 5 cosmological, 6 describing the galaxy–matter connection and its stochasticity, and 2 specifying the intra-halo velocity dispersion. These parameters are listed in Table 1.

For the three main cosmological parameters, we expect to be most constrained by our data |$\lbrace \Omega _m h^2, \sigma _8, h\rbrace$| we will adopt flat priors that coincide with the ranges provided by the emulator. For the two remaining cosmological parameters, |$\lbrace n_s$|⁠, |$\Omega _b\rbrace$| we will adopt informative Gaussian priors given by the marginalized posterior of CMB analyses. For all the remaining (nuisance) parameters, we will adopt flat and uninformative priors. The stochastic parameters |$\epsilon _1$| and |$\epsilon _2$| allow for deviations with respect to Poisson statistics due to overlapping (or contact) effects. Even taking into account that the Poisson noise has been removed from the multipoles, the priors on the stochastic parameters explore non-physical values. We make this choice to absorb any residual observational systematics present in the data at the newly included scales. However, as we will see, we find no compelling evidence for non-physical values.

We note that our bias parameters describe the average relation between galaxies in a given sample and the underlying Lagrangian density field. Therefore, in our fiducial analysis we have an independent set for each BOSS subsample. For instance, when combining the low-z NGC and SGC measurements, we will have 21 free parameters: 5 cosmological and (2|$\times$|8 =)16 independent bias parameters.

It is interesting to note that for a given population, we do not expect the bias parameters to be completely independent. For instance, there are well-established relations among bias parameters for dark matter haloes (see Stücker et al. 2024, and references therein). Although these are not directly applicable to galaxies (due to assembly bias and because galaxy formation physics and selection criteria modulate the haloes entering the sample), the current generation of galaxy formation simulations limits the magnitude of the departures from the halo bias relations. Additionally, as our understanding of galaxy formation improves, these bias relations will become more robust and provide a clear path to incorporate galaxy formation physics to improve cosmological constraints.

To explore this idea, we will complement our fiducial results with an analysis in which we express higher order bias parameters as a function of the linear bias |$b_1$|⁠. Explicitly, we will adopt the mean relationship for Lagrangian bias parameters reported by Zennaro et al. (2022):

which represent the mean over thousands of catalogues of different redshifts, selection criteria, underlying cosmology, and sample number densities. Importantly for consistency, these measurements were obtained by fitting the real-space galaxy power spectrum and the matter-galaxy cross correlation using the same hybrid-bias model of this paper. In Pellejero Ibañez et al. (2022), we checked that these values were consistent between real- and redshift-space.

We will refer to this alternative analysis as ‘Physical Bias Priors’. By adopting information about the relationship among bias parameters, we will reduce the degrees of freedom in our model. This results in a more deterministic prediction for a given cosmological model, lower prior-volume effects, and thus stronger cosmological constraints (see e.g. Ivanov et al. 2024). An study on how strong the choice of such priors need to be to affect the inference of cosmology as performed by Barreira, Lazeyras & Schmidt (2021) is beyond the scope of this paper. By fixing their values, this test shows a ‘best-case’ scenario, where we have independent measurements of the bias coevolution relations. This gives an upper limit to the information that can be extracted with our model using the two-point statistics. Note that further priors could be assumed on other nuisance parameters (see e.g. Schmittfull et al. 2019, 2021; Kokron et al. 2022). However, we make the choice of leaving them uninformative to account for small-scales systematics.

4 VALIDATION ON MOCK CATALOGUES

Before presenting the main results of this study, we validate our methodology by applying it to the Nseries mock catalogues⁵ (Alam et al. 2017).

4.1 The Nseries simulation suite

Nseries is a suite of 84 quasi-independent mock catalogues of the CMASS NGC sample of the BOSS survey. These mocks were generated from seven independent N-body simulations with |$2500\, h^{-1}\,{\rm Mpc}$| a side and adopting the following cosmological parameters: |$\Omega _{\rm m} =0.286$|⁠, |$\sigma _8=0.82$|⁠, |$h= 0.7$|⁠, |$n_s=0.97$|⁠. Mock galaxies were creating by using an HOD formalism, with parameters expected to describe the CMASS sample. Each simulation was projected into 12 different orientations adopting the same radial and angular selection function and fibre incompleteness as the BOSS data (Alam et al. 2017).

Although these mocks only encompass the CMASS NGC sample, spanning |$0.43 \lt z \lt 0.7$| (with an effective redshift of |$z_{\rm eff} = 0.56$|⁠), the substantial volume of the mocks and the detailed modelling of RSDs, allow us to test the ability of our model to deliver unbiased cosmological constraints. Given the differences with respect to the complete BOSS survey, a specific covariance matrix is required for the Nseries. This was measured with 2048 Patchy simulations adopting the window function of the Nseries mocks (Philcox & Ivanov 2022). Finally, the mocks were transformed assuming the fiducial cosmology: |${\Omega _{\rm m,fid}= 0.31,h_{\rm fid}= 0.676}$|⁠, which will test our treatment of geometric RSDs.

4.2 Validation

We validate our model and likelihood by inferring cosmological parameters from the power spectrum multipoles measured in the Nseries mocks.

We use the mean over the 84 mocks as our synthetic data vector. We consider two cases for the uncertainties. First, we employ a covariance matrix corresponding to the total effective volume of BOSS (⁠|$V_{\rm eff} = 6 h^{-3}{\rm Gpc}^3$|⁠). Secondly, we consider a covariance matrix rescaled by a factor of 5, roughly corresponding to the full volume of a DESI-like survey (⁠|$V_{\rm eff} = 30 h^{-3}{\rm Gpc}^3$|⁠). We dub these cases as ‘Nseries’ and ‘DESI’, respectively. The first case will validate the use of our pipeline to BOSS data, whereas the second analysis will provide a much more stringent test of our theoretical modelling.

We display our results in Fig. 1. Our first finding is that we obtain unbiased constraints on the three main cosmological parameters at |$k_{\rm max} =0.4 \, h\, {\rm Mpc}^{-1}$|⁠. In addition, the posterior distributions are consistent when only considering larger scales: the cosmological constraints are statistically compatible for all the values of |$k_{\rm max}$| we consider. Additionally, we observe a progressive increase in constraining power when including small scales, specially in |$\sigma _8$|⁠. For instance, the |$1\sigma$| region on |$\sigma _8$| decreases by |$\sim 25~{{\ \rm per\ cent}}$| when comparing |$k_{\rm max} =0.1$| and |$0.4 \, h\, {\rm Mpc}^{-1}$|⁠. On the other hand, uncertainties on |$\Omega _{\rm m}$| and h are fairly stable for |$k_{\rm max} \gt 0.2\, \, h\, {\rm Mpc}^{-1}$|⁠. As discussed in Pellejero-Ibañez et al. (2020a), we do not expect very large gains beyond |$k_{\rm max} = 0.2 \, h\, {\rm Mpc}^{-1}$| for the volume and number density of tracers in BOSS. Indeed, by using our methodology down to |$0.4 \, h\, {\rm Mpc}^{-1}$|⁠, the figure of merit can increase by a factor of 1.1 compared to the scales usually employed in similar analyses, |$k_{\rm max} =0.2 \, h\, {\rm Mpc}^{-1}$|⁠. This is equivalent to having access to a survey volume 20 per cent times larger.

When considering the ‘DESI’ case (purple contours), we see that our constraints on all parameters agree within the true cosmology (indicated by dashed lines) at the |$1\sigma$| level. We note that these contours coincide with those from the Nseries analysis, which suggests that projection effects in the parameter space are not significant for |$k_{\rm max} =0.4 \, h\, {\rm Mpc}^{-1}$|⁠. However, there can be projection effects at |$k_{\rm max} =0.2 \, h\, {\rm Mpc}^{-1}$|⁠, as indicated by the small shifts in |$\sigma _8$|⁠. Although not shown here, we tested with a DESI volume covariance, and the shifts reduced their statistical significance. We explore this further in Section 6.2.

Therefore, we conclude that theoretical uncertainties in our model are negligible compared to statistical uncertainties in BOSS, and thus our pipeline should deliver robust and accurate parameter constraints using mildly non-linear scales up to |$k_{\rm max} =0.4\, \, h\, {\rm Mpc}^{-1}$|⁠.

5 RESULTS

After having validated our inference pipeline, we now present the key result of this work: parameter constraints from the BOSS DR12 data set. As detailed in previous sections, our analysis relies on the BACCO hybrid Lagrangian bias expansion, which we apply to measurements of the monopole, quadrupole, and hexadecapole of the redshift-space power spectrum in SDSS-III BOSS down to |$0.4\, \, h\, {\rm Mpc}^{-1}$|⁠.

5.1 Cosmological constraints from SDSS-BOSS

In Fig. 2, we present our constraints on the three main cosmological parameters we included in our analysis: |$\Omega _{\rm m}$|⁠, |$\sigma _8$|⁠, and h. We display three cases. The blue and orange contours indicate constraints from the SGC and NGC data, respectively, where we jointly analysed multipoles from the low- and high-z subsamples. The third case, indicated by dashed contours, corresponds to the combination of all four BOSS subsamples modelled adopting physical priors on the relationships of bias parameters (cf. Section 3.2.2).⁶ For comparison, we also display the results obtained from the Planck satellite. The best-fitting parameters as well as the 68 per cent confidence intervals are provided in Table 2. In Appendix A and Fig. A1, we show our constraints on the full parameter space.

Our fiducial NGC analysis yields a 3.6 per cent measurement on |$\Omega _{\rm m}$|⁠, 4 per cent on |$\sigma _8$|⁠, and 2 per cent on h. The accuracy on the derived parameter, |$S_8 \equiv \sigma _8 \sqrt{ \Omega _{\rm m} /0.3}$|⁠, is 4.7 per cent. In the case of the SGC data set, we find significantly larger uncertainties (around 60 per cent larger errors in each of the cosmological parameters), which is expected due to the volume |$\sim 3$| times smaller than NGC.

Notice that the marginalized constraints from NGC and SGC are statistically consistent, with a mild disagreement in the |$\Omega _m$|-h plane. However, there is a better agreement in terms of their best-fitting parameters (shown as filled stars), which suggests the lack of significant sources of unaccounted systematic errors. In turn, the marginalized statistics might be affected (especially for SGC) by the boundaries of our priors and projection effects in the full parameter space.

When considering the case of ‘Physical Bias Priors’ (dashed contours), our estimates show very small shifts compared to those in our fiducial NGC analysis. Furthermore, this case shows a |$\sim 25~{{\ \rm per\ cent}}$| decrease in uncertainties, mostly due to the combination of NGC and SGC.

We have also analysed NGC and SGC separately including our bias priors. The resulting marginalized constraints are shown in grey in Fig. 4. We see that the constraints barely change in the case of NGC, but significantly improve in the case of SGC, due to the reduced freedom in the model.

The consistency in our analyses indicates that the BOSS bias parameters are compatible with our assumed bias relations, as we will show explicitly below. Additionally, the increase in constraining power suggests that bias priors break internal model degeneracies, and are effective in reducing prior volume effects.

Finally, we highlight that both of our models, and therefore the minimal |$\Lambda$|CDM, are a very good description of the BOSS data. We display an example of this in Fig. 3, where we show the multipoles of each BOSS subsample together with the fiducial model evaluated at the best-fitting parameter set obtained from the joint high-z and low-z NGC and SGC analysis. In the bottom panels, we display the difference in units of the diagonal elements of the covariance matrix. In each panel, we provide the value of |$\chi ^2_{\rm red} = \chi ^2/{\rm d.o.f}$|⁠, where for simplicity we have estimated the number of degrees of freedom |${\rm d.o.f}$| as the number of data points minus the number of free parameters. The values are very consistent with 1, showing the ability of the model to fit the range of scales considered in this work. Although not shown, the |$\chi _{\rm red}^2$| value increases slightly for the ‘physical bias prior’ case by approximately |$\Delta \chi _{\rm red}^2 \sim 0.05$|⁠.

5.2 Comparison to previous BOSS analyses

The BOSS data set is the largest publicly available galaxy redshift survey, and has been previously analysed by a large number of works (see d’Amico et al. 2020; Ivanov et al. 2020; Tröster et al. 2020 for some of the first examples on this ‘full-shape’ analysis). Fig. 4 provides a compilation of recent constraints on |$\Omega _{\rm m}$|⁠, |$\sigma _8$|⁠, h, and |$S_8$| using similar statistics to this work.

More specifically, the same datavector we employ here was analysed by Ivanov et al. (2020) and Philcox & Ivanov (2022). They used EFT to model the power-spectrum multipoles down to |$k_{\rm max} \simeq 0.2\, \, h\, {\rm Mpc}^{-1}$|⁠. D’Amico et al. (2024) used a similar model, but on the windowed BOSS data (Beutler & McDonald 2021) up to |$k_{\rm max} =0.23\, \, h\, {\rm Mpc}^{-1}$|⁠. Finally, Kobayashi et al. (2022) analysed the data with an HOD emulator up to |$k_{\rm max} = 0.25 \, h\, {\rm Mpc}^{-1}$|⁠, also with the windowed version of the datavector. As in our analyses, these works adopt informative priors on |$n_s$| based on Planck and on |$\Omega _bh^2$| based on BBN (Cooke, Pettini & Steidel 2018).

Overall, we find that our constraints agree statistically with those obtained in previous works. However, our analysis delivers |$\sigma _8$| measurements which are 30 per cent and 40 per cent more accurate in the fiducial NGC and physical bias prior cases, respectively. There is, however, considerably scatter among the measurements of |$\sigma _8$|⁠, as seen in the second column of Fig. 4. This could suggest that the specific choices of each analysis plays an important role.

Specifically, previous works based on EFT of LSS (Carrasco et al. 2014) used |$k_{\rm max} \simeq 0.2\, \, h\, {\rm Mpc}^{-1}$|⁠, since the model accuracy drops rapidly on smaller scales. In contrast, our model has been shown to work well on mocks down to scales of |$k\simeq 0.6 \, h\, {\rm Mpc}^{-1}$|⁠. There are two possible differences between EFT approaches and our approach. On one side, we are pushing our study to smaller scales which potentially carry information. As we will see in Section 6.1, the increase on constraining power with scale in our fiducial case is negligible, therefore making this explanation highly unlikely. On the other side, our model does not need the use of ‘counterterms’, since all the non-linear evolution is encoded in the N-body displacement field.⁷ This changes the shape of the parameter space, changing the degeneracies between cosmological and nuisance parameters.

The derived structure parameter |$S_8$| is also consistent among different analyses. The biggest disagreement is with D’Amico et al. (2024), who reported a systematically higher value for |$S_8$| than every other BOSS analysis. This result, however, includes a correction for prior volume effects, which makes it difficult to compare against other measurements. In contrast, our constraints at |$k_{\rm max} \simeq 0.4 \, h\, {\rm Mpc}^{-1}$| do not show evidence of strong prior volume effects, as we will discuss in Section 6.2.

The BOSS multipoles were analysed jointly with additional summary statistics by Philcox & Ivanov (2022), Chen et al. (2022b), Ivanov et al. (2023), and D’Amico et al. (2024), who included combinations of the galaxy bispectrum monopole, quadrupole, an analogue to the real-space power spectrum, and the post-reconstruction BAO feature. Consistent with these works, our results are in agreement in terms of the marginalized statistics and remain among the most precise determinations of |$\sigma _8$| despite using a reduced data vector. For instance, the measurements of Ivanov et al. (2023), which incorporated all of the aforementioned statistics, yield |$H_0=68.2\pm 0.8$| km s|$^{-1}$|Mpc|$^{-1}$|⁠, |$\Omega _{\rm m} = 0.33\pm 0.01$|⁠, |$\sigma _8=0.736\pm 0.033$|⁠, and |$S_8= 0.77\pm 0.04$|⁠. Compared to our results, their analysis shows 38 per cent stronger constraints for |$H_0$|⁠, arguably due to the post-reconstruction BAO, but slightly looser constraints for |$\Omega _{\rm m}$| and |$\sigma _8$| (14 per cent and 24 per cent, respectively).

Another interesting comparison is with the results of Kobayashi et al. (2022). Using an emulator for the redshift-space power spectrum of haloes combined with the HOD model (Kobayashi et al. 2020), they were able to analyse scales up to |$k_{\rm max} \approx 0.3 \, h\, {\rm Mpc}^{-1}$|⁠. Although our approach to include RSD in the model differs, it is similarly based on halo velocities. Despite differences in our modelling strategies, we find consistent results. Notably, our method yields a slight increase in accuracy for |$\sigma _8$|⁠, achieving 4 per cent accuracy compared to their 5 per cent. This is a minor difference given the significantly different analysis methods. None the less, both studies highlight the potential of using N-body emulator-based models to extract cosmological information.

Finally, it is noteworthy that our constraints on h and |$\Omega _{\rm m}$| align remarkably well with those obtained by the latest BAO results from the DESI collaboration. The values presented in DESI Collaboration (2024) are |$\Omega _{\rm m}=0.295 \pm 0.015$| for DESI alone and |$h=0.6929 \pm 0.0087$| when combined with |$r_d$| from the CMB. This consistency follows the trend towards lower values of |$\Omega _{\rm m}$| and higher values of h observed in Fig. 2.

5.3 Constraints on galaxy bias parameters

The Lagrangian bias parameters measure how, on average, the number of galaxies is modified by large-scale features of the Lagrangian density field. For instance, |$b_1$| measures how galaxies respond to large-scale changes in density, whereas |$b_{s^2}$| does so for changes in the tidal field. In this sense, the bias parameters have a clear physical meaning, and it is therefore interesting to compare their values in observations to theoretical expectations, as this is a direct test of galaxy formation models.

In Fig. 5, we show our constraints on the Lagrangian bias parameters from the analysis of the redshift-space multipoles, using |$k_{\rm max} =0.4 \, h\, {\rm Mpc}^{-1}$|⁠. Note that we display our results in terms of the linear bias, |$b_1$|⁠, since the absolute value of the bias parameters is much more sensitive to cosmology and details of the galaxy population than their relation to |$b_1$| is. We display results from the analysis of BOSS data and of the Nseries mocks. For comparison, we also display the relationship measured by Zennaro et al. (2022) from a large suite SHAMe galaxy mocks, exploring widely different galaxy formation parameters, redshifts, cosmologies, and number densities, including one sample that resembles the clustering of galaxies stellar-mass selected galaxies in the IllustrisTNG simulation with the same number density as BOSS CMASS sample.

We find that the recovered bias parameters agree with the Zennaro et al. (2022) relationship within the expected uncertainty. In particular, the Nseries sample is within |$0.4\sigma$| of the mean predicted value, albeit the |$b_1$| value being systematically smaller than the BOSS high-z data (this is in part because the Nseries simulations adopted a 15 per cent higher |$\sigma _8$| than our BOSS best-fitting value). The bias parameters in BOSS show more scatter than in Nseries (since for the latter we employ the mean of 84 mocks) but always lie within the predicted expected region and very close to the mean relation (shown by solid black lines).

The fact that the prediction roughly agrees with the measured values indicates that, even though the bias parameters might absorb some of the observational systematic not taken into account in the observational weights, they seem to retain their physical meaning as response functions to the large-scale density field.

The only subsample that shows somewhat of a tension with the theoretical expectations is the SGC high-z sample. In fact, this sample not only deviates from the mean |$b_1{\!-\!}b_2$| and |$b_1{\!-\!}b_{s^2}$| relationships, but the value of |$b_1$| is also different from that of the NGC high-z despite both samples containing similar types of galaxies. Although this tension is not strong from a statistical point of view, we note that when we adopt bias priors in the SGC analysis, the derived value is |$b_1=1.49$|⁠, in much better agreement with the NGC counterpart.

6 ROBUSTNESS TESTS

In this section, we explore the robustness of our results to various analysis choices. Specifically, we quantify the sensitivity of our cosmological constraints to |$k_{\rm max}$|⁠, the importance of the hexadecapole, and the impact of prior volume effects.

6.1 The impact of |$k_{\rm max}$|

In Fig. 6, we compare the marginalized constraints on |$\Omega _{\rm m}$|⁠, h, and |$\sigma _8$| as a function of the maximum scale included in the analysis, |$k_{\rm max}$|⁠. We display our fiducial analysis carried out independently on NGC and SGC data, and of their joint analysis in our ‘physical bias priors’ case.

First, we see that our constraints remain remarkably stable from |$k_{\rm max} =0.2$| to |$0.4\, \, h\, {\rm Mpc}^{-1}$| (note we do not consider |$k_{\rm max} \le 0.1\, \, h\, {\rm Mpc}^{-1}$| as these constraints are heavily affected by the limits adopted in our emulator). In almost all cases, the parameter shifts are below |$0.4\sigma$|⁠, with the exception of h in the NGC case, which shifts by |$0.7\sigma$|⁠. However, h is extremely stable when we adopt bias priors, which suggest that, at least part of this can be explained by internal model degeneracies at |$k_{\rm max} =0.2\, \, h\, {\rm Mpc}^{-1}$|⁠. This is also evident from the comparison of the best-fitting values, shown as star symbols, which display a higher level of consistency than the marginalized statistics in the fiducial case. This indicates that the best fit obtained at |$k_{\rm max} =0.4, \, h\, {\rm Mpc}^{-1}$| also provides a good fit at |$0.2, \, h\, {\rm Mpc}^{-1}$|⁠. We verified this and found a value of |$\chi ^2_{\rm NGC,red}(k_{\rm max} =0.2, \, h\, {\rm Mpc}^{-1}) \approx 1.04$|⁠.

Secondly, we note that the accuracy on cosmological parameters does not improve significantly with |$k_{\rm max}$| in our fiducial NGC and SGC analyses. This is consistent with our results on the Nseries mocks (cf. Section 4) and with the findings of Pellejero Ibañez et al. (2024), where we used the BACCO hybrid bias emulator to analyse a mock catalogue resembling BOSS. Due to the limited volume and number density of BOSS galaxies, the additional small-scale Fourier modes mostly constrain the nuisance parameters of the model. For instance, the accuracy with which the bias parameters are measured increases by a factor of 30 per cent from |$k_{\rm max} =0.2$| to 0.4 (mainly in |$b_2$| and |$b_{s^2}$|⁠), and by 90 per cent in the case of the parameters modelling small-scale velocities (⁠|$\lambda _{\rm FoG}$| and |$f_{\rm sat}$|⁠).

Given that small scales in BOSS mostly constrain nuisance parameters of our model, it is not surprising that in the ‘physical bias priors’ case we detect a mild but steady increment in cosmological information with |$k_{\rm max}$|⁠. For instance, the estimate of |$\sigma _8$| improves by 32 per cent compared to |$k_{\rm max} =0.2 \, h\, {\rm Mpc}^{-1}$|⁠. However, the gains are still moderate considering the enormous increase in the number of modes. This is because of the additional model freedom provided by the nuisance parameters describing RSDs and the stochasticity in the matter−halo connection. This highlights the potential gains of investigating and placing robust priors on these parameters, in an analogous way to the bias parameters.

6.2 Prior volume effects

Due to the complexity of our model, there are several nuisance parameters that could be loosely constrained, especially when considering low values of |$k_{\rm max}$|⁠. In such cases, the marginalized constraints on cosmological parameters might suffer from ‘prior-volume effects’, that is, the estimated posterior becomes sensitive to the prior distribution, boundaries adopted, and projections from a high-dimensional parameter space. In the context of EFT of LSS, Carrilho, Moretti & Pourtsidou (2023) and D’Amico et al. (2024) have cautioned that the analysis of BOSS might suffer from such projection effects. Although the relevance of priors is a key aspect of Bayesian statistics, it is important to identify the existence of such effects, otherwise it might mislead the interpretation of results.

In Fig. 7, we present a quantification of prior-volume effects in our results. Specifically, we present the analysis of a datavector generated by evaluating our model with the best-fitting values of the combined BOSS-NGC sample (see Table 2). The data vector is then noiseless and represents the mean of all possible realizations of the theory as defined by the covariance from the Patchy mocks. Naturally, we expect the best-fitting parameters to exactly coincide with the input values (indicated by dotted lines). Therefore, any deviation between the input and the marginalized values can be attributed to the effect of priors.

We can see that projection effects are almost negligible in our fiducial analysis of NGC, implying our current priors are uninformative. This also eliminates the need for correcting our results in a manner explored by D’Amico et al. (2024). In contrast, when the model is used only up to |$k_{\rm max} =0.2~ \, h\, {\rm Mpc}^{-1}$|⁠, the mean of the marginalized posterior on |$\sigma _8$| and |$\Omega _{\rm m}$| overestimate the true values by almost |$1\sigma$|⁠, as a consequence of poorly constrained nuisance parameters. Note that this does not entirely demonstrate that we are unaffected by these effects; rather, it shows that we are not dominated by them at the maximum likelihood values of our parameter space. Different regions will be affected differently by prior-volume effects. To completely rule out being dominated by them, a profile likelihood analysis would be required. However, this is computationally expensive, and we leave such a test for future work.

6.3 The impact of the hexadecapole

The hexadecapole of the galaxy power spectrum is the most sensitive to the line-of-sight correlations, and hence to the physics causing small-scale velocities. In addition, it is the noisiest multipole predicted by our emulator, owing to the finite volume of the simulations employed. Additionally, in Beyond-2pt Collaboration (2024), we showed that the uncertainties in the emulation of the hexadecapole could introduce small but significant biases in our estimation of |$\sigma _8$|⁠.

To explore the robustness of our results to the modelling of the hexadecapole, we re-analysed the BOSS-NGC data, but only considering the monopole and quadrupole. In Fig. 8, we present our results for |$k_{\rm max} =0.2$| and |$0.4\, \, h\, {\rm Mpc}^{-1}$|⁠. In both cases, we see that the hexadecapole shifts the inferred parameter values by less than |$0.2\sigma$| – a change that is statistically insignificant compared to other sources of noise. This is consistent with our expectations from the analysis of the Nseries mocks, which showed no significant parameter biases, and with Pellejero Ibañez et al. (2023) where we showed that the emulator noise is below the statistical uncertainty of BOSS data.

Finally, it is worth noting that although the hexadecapole has a negligible impact on the precision when |$k_{\rm max} =0.4\, \, h\, {\rm Mpc}^{-1}$|⁠, it improves the accuracy in the |$k_{\rm max} =0.2 \, h\, {\rm Mpc}^{-1}$| case, where, for instance, the uncertainty in |$\sigma _8$| decreases by 10 per cent. Additionally, when excluding the hexadecapole, the cosmological constraints consistently improve with |$k_{\rm max}$|⁠, unlike in our fiducial analysis (cf. Section 6.1). This indicates that at least part of the information encoded on small scales is already captured by the hexadecapole on large scales.

7 DISCUSSION

7.1 The |$S_8$| tension

It has become increasingly common to characterize the amplitude of matter fluctuations in the Universe via the ‘structure parameter’ |$S_8$|⁠. The motivation behind this is that weak gravitational lensing directly measures this parameter combination, which has influenced its use in other low-z cosmic probes. Comparing direct measurements of |$S_8$| at low redshift, with the expectations of high-z measurements offers a direct way to test whether the growth of structure is compatible with the predictions of the |$\Lambda$|CDM model.

With the arrival of Stage-III lensing surveys, a tension started to emerge between low-z measurements of |$S_8$| and the expectation from analyses of the Planck satellite. For instance, the analysis of shear correlations as measured by DES-Y3, HSC-Y3, and KIDS1000, all reported a tension with Planck, ranging between 2 and |$3\sigma$| (Asgari et al. 2021; Heymans et al. 2021; Secco et al. 2022; Dalal et al. 2023). Similarly, the joint analysis of lensing with photometric galaxy clustering (a.k.a. 3|$\times$|2pt analyses), strengthened this conclusion, typically reporting a slightly larger tension (Heymans et al. 2021; Abbott et al. 2022). One of the most precise measurements of |$S_8$| was provided by García-García et al. (2021), who combined multiple surveys to find a value |$3.4\sigma$| smaller than in Planck. Similarly, analysis of the cross-correlation between galaxies and CMB lensing, and of the abundance of clusters, as well as the analyses of BOSS clustering, all reported values systematically below Planck (Bocquet et al. 2019; Ivanov et al. 2020; Krolewski, Ferraro & White 2021; White et al. 2022). Additionally, Amon et al. (2023) suggested that a low |$S_8$| value could be the solution for the ‘lensing-is-low’ problem (Leauthaud et al. 2017). The state was summarized in a 2021 Snowmass report (Di Valentino et al. 2021).

More recently, the situation has started to change. As emphasized in an independent reanalysis of DES-Y3, Aricò et al. (2023b) cautioned that current weak lensing constraints were sensitive to the model for the non-linear power spectrum, the treatment of intrinsic alignments, and the assumptions regarding baryonic physics. Explicitly, using the BACCO emulators (Aricò et al. 2021, 2023a), Aricò et al. (2023b) found a value for |$S_8$| compatible with Planck at the |$1.4\sigma$| level. This was also the conclusion of a subsequent official reanalysis of the DES and KIDS surveys (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023). Additionally, an important point to note is that most of the constraining power of weak lensing arises from non-linear scales. Therefore, as highlighted by Amon & Efstathiou (2022), agreement between lensing and Planck can be obtained if baryonic physics is stronger than the level that is typically assumed in hydrodynamical simulations.

An alternative way of measuring gravitational lensing is through its effect on the CMB: this signal is well consistent with the prediction of the Planck cosmology (Planck Collaboration VI 2020; Madhavacheril et al. 2024). However, the total CMB lensing signal comes from much higher redshifts than those probed by weak galaxy shear. A tomographic study that isolated the |$z\lt 0.8$| part of the CMB lensing signal found a result lower than Planck (Hang et al. 2021), arguing that this could be consistent with both the low |$S_8$| results and the good agreement between total CMB lensing and Planck if the tension was mainly in the direction of lower |$\Omega _m$|⁠.

Finally, the most recent analyses of the abundance of clusters detected by the SPT survey (Bocquet et al. 2024) and of the cross correlation of Planck-CMB lensing with galaxies in the WISE survey (Farren et al. 2024) both reported milder tensions with Planck. The eROSITA cluster abundances (Ghirardini et al. 2024) even yielded values of |$S_8$| slightly above that of Planck, with a nominal accuracy that even exceeds that of Planck. Finally, Chaves-Montero, Angulo & Contreras (2023), Contreras et al. (2023a), and Contreras, Chaves-Montero & Angulo (2023b) showed that an apparent ‘lensing-is-low’ tension appears due to the limitations of simplistic models for the galaxy–halo connection.

In Fig. 9, we show a compilation of recent constraints on |$S_8$|⁠, compared with our BOSS results (shown in the first three rows). We can see that our results are statistically compatible with most of the recent low-z cosmic probes, perhaps with the exception of eROSITA. Moreover, our reported |$S_8$| value is |$2{\!-\!}3\sigma$| lower than Planck with a mean value towards the lower end of values in the literature, although still compatible with other full-shape analyses such as Ivanov et al. (2023). Overall, the |$S_8$| tension is statistically mild: given the scatter between the non-Planck measurements, it is hard to make an unambiguous case that there is an inconsistency that requires systematics. Nevertheless, almost all alternative measurements argue that |$S_8$| should be below the central Planck value, and so at a minimum there is a good case that the Planck value has fluctuated high.

In the next subsections, we will explore the |$S_8$| tension further, first in terms on an alternative parametrization for the amplitude of structure fluctuations, |$S_{12}$|⁠, and then on the feasibility of the Planck cosmology to explain BOSS data.

7.2 |$S_{12}$| and the |$H_0$| tension

In light of the tension between our best-fitting parameters with Planck, in this section we explore a different parametrization and the role of an informative prior on |$H_0$|⁠.

Motivated by Sánchez (2020), we consider the parameters |$\omega _m \equiv \Omega _m h^2$| and |$S_{12} = \sigma _{12} (\omega _m/0.14)^{0.4}$|⁠, where |$\sigma _{12}$| is the rms linear fluctuations in 12 Mpc spheres. The advantages of this parametrization are that the value of h does not enter in the definition of |$S_{12}$| and that there is a separation into parameters that affect the shape and those that affect the amplitude of the linear power spectrum. Furthermore, as discussed in Sánchez (2020) and García-García et al. (2024), weak lensing constrains the overall amplitude of fluctuations which makes it sensitive to the adopted prior on h. In fact, García-García et al. (2024) showed that the tension in |$S_8$| between weak lensing and Planck does not appear either in |$S_{12}$| or when a prior on |$H_0$| based on SH0ES data is employed, suggesting that the current |$S_8$| tension might be another face of the |$H_0$| tension between Planck and SH0Es.

In Fig. 10, we display our constraints on |$\omega _m$| and |$S_{12}$| from our fiducial NGC analysis. Unlike in weak lensing analyses, adopting this parametrization does not relax the tension between BOSS and Planck. Additionally, imposing a prior on |$H_0$| from local measurements from SH0ES (⁠|$h = 0.7304 \pm 0.0104$|⁠, Riess et al. 2022) does not alleviate the tension with Planck – in fact, the discrepancy in |$S_{12}$| increases slightly. The main reason behind these results is that, unlike weak lensing analyses, full-shape clustering already constrains |$H_0$| very well (⁠|$\sim 2~{{\ \rm per\ cent}}$|⁠, similar to the accuracy of SH0ES), which implies an almost unique relationship between |$S_{12}$| and |$S_8$|⁠. Hence, a prior on |$H_0$| will only have a moderate effect on the constraints from redshift-space clustering.

7.3 Comparison with Planck cosmology

In Section 3.2, we showed |$\Lambda$|CDM to be an excellent description of the data, however, as can be seen in Figs 2 and 9, the parameter constraints differ with those obtained from the Planck satellite.

We quantify the tension between BOSS and Planck in the full cosmological space as

where |$\mu \equiv \lbrace \Omega _m, \sigma _8, h\rbrace$| and C is the marginalized parameter covariance, and the subscripts BOSS and |${\it Planck}$| refer to quantities estimated by our BOSS analysis and by Planck, respectively. In this three-dimensional space, we find a 2.3|$\sigma$| tension between Planck and BOSS, which remains at the same level in our ‘physical bias prior’ analysis (a decrease in the tension in h owing to SGC data is compensated by increased tensions in |$\Omega _{\rm m}$| and |$\sigma _8$|⁠).

To explore this tension further, we reanalysed the BOSS data but this time fixing the cosmological parameters to the values preferred by Planck. This will tell us whether there exists a set of bias and nuisance parameters that provide a good fit to BOSS within the Planck cosmology.

Indeed, in Fig. 11 we see that Planck and our best-fitting cosmology both provide a good description of the low and high-z BOSS data, especially at |$k \lt 0.25\, \, h\, {\rm Mpc}^{-1}$|⁠. Although not shown here, the two parameters showing the greatest deviations when fixing to this new cosmology are |$b_1$| and |$b_{s^2}$| (which are also more tightly constrained). Specifically, the linear bias |$b_1$| decreases to accommodate for a higher value of |$\sigma _8$|⁠, while the tidal field bias parameter shows a preference towards positive non-zero values, |$b_{s^2}\simeq 1.3$| at the 2|$\sigma$| confidence level. Note that from our studies in Zennaro et al. (2022), it is not possible to completely rule out the possibility for BOSS galaxies to live in regions where the Lagrangian tidal field is positive, specially given the complex selection criteria of BOSS and that these bias relations were built for a generic galaxy formation sample.

The value of |$\chi ^2$| increases only mildly when the cosmology is fixed to Planck, from 446.03 to 448.02. From a frequentist point of view, the difference of |$\Delta \chi ^2=2$| is expected to occur 60 per cent of the times considering that we have 456 data points, minus some allowance for the non-linear degrees of freedom in the model, which we neglect. Thus, apparently the Planck cosmology is in good agreement with the BOSS data when looked at in this way. This may however partly reflect the fact that the |$\chi ^2$| test is something of a blunt tool: it lacks knowledge of the form of the model under study, which is an explicit feature of the Bayesian analysis.

Regarding the degree to which the model’s non-cosmological parameters need to be fine-tuned to accommodate Planck, we can assess this by computing the Bayesian evidence ratio, R, defined as the ratio between the evidences of a model adopting the Planck set of values and our BOSS cosmology. The result is |$\ln R \simeq -0.5$|⁠, which is well below the commonly accepted threshold of 4.6 for a decisive preference for one model over another (Jeffreys 1998). Hence, we conclude that the full-shape BOSS clustering is in no significant tension with the Planck set of cosmological parameters from the Bayesian point of view, despite the decrease in number of dimensions.

Finally, as shown in dotted lines in Fig. 11, we can repeat this comparison using our ‘physical bias prior’ analysis. By using the Planck cosmological parameters and further fixing the high-order bias parameters to their ‘physical bias prior’ values, we find a 1|$\sigma$| distance from the fiducial result in the best-fitting models, |$\Delta \chi ^2 = 15$|⁠, and a Bayes ratio of |$\ln R \simeq -8.2$|⁠, indicating a decisive evidence in favour of this case over that where cosmology is varied, despite the larger |$\chi ^2$|⁠. This result can be interpreted as follows: the nuisance parameters of our model can accommodate a fit where the cosmological parameters are set to Planck and where the ‘physical bias prior’ are used. This fit is only partially excluded, since it is only 1|$\sigma$| away from the fiducial result. However, since the dimensionality reduced from 21 to 10, the Bayes factor heavily penalizes the inclusion of these extra free parameters.

A fairer test is to compare cases where the cosmology is either fixed to Planck or to the best-fitting BOSS-NGC cosmology. This gives a Bayes ratio of |$\ln R \simeq 2.1$|⁠, providing substantial evidence for the BOSS cosmology over Planck, when assuming ‘Physical Bias Priors’. We caution, however, that we do not expect BOSS galaxies to exactly coincide with our bias priors. Instead, this result highlights the benefits of adopting informative bias priors: although our cosmological constraints did not improve significantly after adopting relationships among bias parameters, such priors can be extremely valuable when comparing competing models. Therefore, we emphasise the future need for robust priors on nuisance parameters, estimated according to our state-of-the-art understanding of galaxy formation physics.

8 CONCLUSIONS

In this paper, we have presented a novel analysis of the multipoles of the redshift-space power spectrum of galaxies in the SDSS-IV BOSS survey. We employed a novel theoretical model based on a perturbative description of the relationship of matter and galaxies, which we combined with displacement and velocity fields extracted from cosmological simulations. The advantages of our novel modelling allow us to employ much smaller scales (⁠|$k_{\rm max} \simeq 0.4 \, h\, {\rm Mpc}^{-1}$|⁠) than those typically used previous analyses.

Our fiducial analysis of the NGC provides one of the most accurate measurements to date on |$\Omega _{\rm m}$|⁠, h, and |$\sigma _{8}$| from the BOSS data set multipoles alone. Specifically, we obtain an accuracy of 3.6 per cent on |$\Omega _m=0.301\pm 0.011$|⁠, 4 per cent on |$\sigma _8=0.745^{+0.028}_{-0.035}$|⁠, 2 per cent on |$h=0.705\pm 0.015$|⁠, and 4 per cent accuracy on the structure parameter, |$S_8 = 0.747^{+0.032}_{-0.039}$|⁠.

Our constraint on |$S_8$| is |$2.5\sigma$| lower than the value inferred by the analysis of the Planck satellite data. This low |$S_8$| value is compatible with constraints from other low-z probes such as recent analyses of the combination of gravitational shear, angular clustering, and their cross-correlation. However, we argue that the current BOSS data are not sufficient to confidently identify a tension with Planck. Specifically, we showed that the Planck cosmology is still able to achieve a good fit to the power spectrum data, and the differences with respect to our best-fitting cosmology can be absorbed by the astrophysical nuisance parameters of our model.

Another novel aspect of our analysis was the adoption of physical priors on the relationships between the linear and higher order bias parameters. Adopting these priors reduces prior volume effects and increases the accuracy of our fiducial cosmological constraint on |$S_8$| by 10 per cent. However, we note that our priors correspond to the mean expected in a generic galaxy sample, and were not built specifically for a BOSS-like sample. In this sense, this analysis can be seen as exploring the potential gains and benefits of combining flexible models with our information arising from our current understanding of galaxy formation.

In the future, by considering only specific astrophysical scenarios appropriate to the particular galaxy sample under study, we expect to obtain more robust and even more informative priors on the nuisance parameters of our model. Importantly, this approach provides a clear and systematic route to incorporate state-of-the-art galaxy formation simulations into cosmological inference. We plan to explore this in forthcoming publications.

ACKNOWLEDGEMENTS

MPI is supported by the Science and Technology Facilities Council (STFC) consolidated grant no. RA5496. REA acknowledges support from project PID2021-128338NB-I00 from the Spanish Ministry of Science and support from the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement number 101086388 (LACEGAL), and the support of the European Research Council Starting Grant (ERC-StG) number 716151 (BACCO). We would like to thank Oliver Philcox for providing with the Patchy-Nseries data, Naomi Robertson, Martin White, Giovanni Arico, David Alonso, Matteo Zennaro, Alex Amon, Masahiro Takada, and Pedro Gregorio Carrilho for useful discussions. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any author accepted manuscript version arising from this submission.

DATA AVAILABILITY

This work used data of the BOSS-SDSS collaboration,⁸ specifically, in its window-free power spectra version.⁹ It also made extensive use of the Nseries mock catalogues.¹⁰ The main theory package used in this work is the baccoemu package,¹¹ and for the MCMC chains, the multinest¹² code.

Footnotes

REFERENCES

APPENDIX A: APPENDIX A

In this appendix, we present the full marginalized posteriors for all cosmological and non-cosmological parameters sampled within the MCMC likelihoods, using the power spectrum multipoles. We display them for the three cases shown in Fig. 2, namely using the BOSS-NGC, BOSS-SGC, and combined BOSS data sets with physical bias priors.