Cosmological constraints from density-split clustering in the BOSS CMASS galaxy sample

We present a clustering analysis of the BOSS DR12 CMASS galaxy sample, combining measurements of the galaxy two-point correlation function and density-split clustering down to a scale of $1\,h^{-1}{\rm Mpc}$. Our theoretical framework is based on emulators trained on high-fidelity mock galaxy catalogues that forward model the cosmological dependence of the clustering statistics within an extended-$\Lambda$CDM framework, including redshift-space and Alcock-Paczynski distortions. Our base-$\Lambda$CDM analysis finds $\omega_{\rm cdm} = 0.1201\pm 0.0022$, $\sigma_8 = 0.792\pm 0.034$, and $n_s = 0.970\pm 0.018$, corresponding to $f\sigma_8 = 0.462\pm 0.020$ at $z \approx 0.525$, which is in agreement with Planck 2018 predictions and various clustering studies in the literature. We test single-parameter extensions to base-$\Lambda$CDM, varying the running of the spectral index, the dark energy equation of state, and the density of massless relic neutrinos, finding no compelling evidence for deviations from the base model. We model the galaxy-halo connection using a halo occupation distribution framework, finding signatures of environment-based assembly bias in the data. We validate our pipeline against mock catalogues that match the clustering and selection properties of CMASS, showing that we can recover unbiased cosmological constraints even with a volume 84 times larger than the one used in this study.


INTRODUCTION
Within the vast cosmic structures we observe today, signatures of primordial features are intermixed with non-linear processes that shape the evolution of galaxies, constituting a ground that is challenging to model but rich in information.In our standard cosmological model, the ΛCDM paradigm, the present-day Universe is the result of a hierarchical structure formation scenario that started from primordial ★ E-mail: enrique.paillas@uwaterloo.cadensity perturbations, which were amplified during inflation and continued to grow through gravitational collapse until today (Guth & Pi 1982;Hawking 1982).Galaxy clustering has been pivotal in testing this hypothesis from late-time Universe data by characterizing the way in which matter is distributed in space, using galaxies as biased tracers of the underlying dark matter distribution.
The most common way to approach galaxy clustering is through the two-point correlation function, or its Fourier pair, the power spectrum.These statistics provide a nearly-complete description of the galaxy density field on large scales and encode information about physics of the early Universe (Peebles 1980).Baryon acoustic oscillations (BAO), which originate from sound waves in the photon-baryon plasma before recombination (Peebles & Yu 1970;Sunyaev & Zeldovich 1970), leave an imprint on the matter distribution, which is detected as a bump in the correlation function or as wiggles in the power spectrum (Eisenstein & Hu 1998).This acoustic feature works as a standard ruler, allowing us to measure the expansion rate of the Universe at different epochs (Percival et al. 2001;Eisenstein et al. 2005;Cole et al. 2005;eBOSS Collaboration et al. 2020;Moon et al. 2023).In addition to BAO, the full shape of the galaxy correlation function and power spectrum contains an enormous wealth of information due to its sensitivity to the growth rate of cosmic structure (Blake et al. 2011;Reid et al. 2012;Alam et al. 2017;Brieden et al. 2021), dark energy, the physics of neutrinos (Zhang & Cai 2022), primordial non-Gaussianity (Moradinezhad Dizgah et al. 2021), and the galaxy-halo connection (Yuan et al. 2022a).
Models based on perturbation theory provide an accurate description of galaxy clustering data on linear and mildly non-linear scales, allowing the extraction of information from the full shape of the power spectrum (e.g., Sánchez et al. 2017;Grieb et al. 2017;d'Amico et al. 2020;Tröster et al. 2020;Bautista et al. 2021;Philcox & Ivanov 2022;Semenaite et al. 2022Semenaite et al. , 2023)).The assumptions behind these models tend to break down in the highly non-linear regime (≲ 20 ℎ −1 Mpc), which has motivated the development of models calibrated on N-body simulations that can accurately describe the galaxy field on scales where non-linear physical processes become relevant.This has recovered information from a portion of the survey data that is usually discarded from the standard clustering analyses, providing clues not only about cosmology but also in regards to galaxy evolution and its connection to the dark matter halo field (Kobayashi et al. 2020;Zhai et al. 2023;Chapman et al. 2022;Lange et al. 2022Lange et al. , 2023;;Yuan et al. 2022a).
Two-point functions provide a complete description of Gaussian density fields.This is satisfied at large scales, where non-linear evolution is mild, and the Gaussianity of the primordial fluctuations are still largely preserved.The late-time Universe, however, is highly non-Gaussian at small scales due to non-linear evolution, and higherorder summary statistics are required to capture all the information.Thanks to substantial theoretical and algorithmic development over the last years, measurements of N-point correlation functions (Slepian et al. 2017;Philcox et al. 2021;Sugiyama et al. 2023) and polyspectra (Gil-Marín et al. 2017;Gualdi et al. 2021;Philcox & Ivanov 2022) are now being performed on data, which not only tightens the constraining power on ΛCDM, but also opens an avenue for testing potential signatures of physics beyond our fiducial model, such as parity violation (Philcox 2022;Hou et al. 2023a).However, the measurement of these N-point statistics remains challenging due to their high computational demands and the large volumes that are needed to detect them with enough statistical significance.This has motivated the development of robust, informative, and efficient clustering methods that can access the information that leaks into higher orders, which can then be complemented and cross-validated with the standard N-point clustering analysis.
Several alternative summary statistics that meet these criteria have been proposed in the literature, including k-th nearest neighbor statistics (Banerjee & Abel 2021), wavelet scattering transforms (Valogiannis & Dvorkin 2022a,b), void statistics (Lavaux & Wandelt 2012), marked correlations (White 2016;Massara et al. 2022), skew spectra (Hou et al. 2023b), and Minkowski functionals (Lippich & Sánchez 2021).Among these novel methods, Paillas et al. (2021) proposed to perform a clustering analysis split by local density, combining the information content of different environments of the cos-mic web.Paillas et al. (2021) showed that density-split clustering can tighten the constraints on geometry and growth by modelling redshift-space distortions around different environments, compared to the standard two-point clustering analysis.Paillas et al. (2023) expanded on this by quantifying the information content of the full shape of the density-split correlation functions, forecasting that the method can deliver precise constraints on ΛCDM parameters, and can potentially be used to put upper limits on the sum of neutrino masses.
Until now, a model that can capture the cosmological dependence of the full shape of the density-split correlation functions was not available.In Cuesta-Lazaro et al. (2023), we have presented sunbird, a simulation-based model for density-split and two-point clustering that can operate down to intra-halo scales, well into the non-linear regime, which has been validated against high-fidelity mock galaxy catalogues based on the AbacusSummit suite of simulations.In this work, we use sunbird to carry out the first application of density-split clustering to observational data, applying it to the final data release of the Baryon Oscillation Spectroscopic Survey (Dawson et al. 2013).We fit a full-shape model of the galaxy two-point correlation function and density-split multipoles down to 1 ℎ −1 Mpc, putting constraints on the base-ΛCDM model, as well as on extensions to the base model that vary the dark energy equation of state, the density of relic neutrinos, and the running of the spectral index of the primordial power spectrum.
The paper is organized as follows.We define our observables in Sect. 2. The clustering modelling, as well as the simulations used to validate our pipeline are presented in Sect.3. We present our main cosmological constraints and the validation tests in Sect. 4. Finally, we summarize and conclude in Sect. 5.

BOSS CMASS galaxy sample
We use data from the final data release (DR12) of the Baryon Oscillation Spectroscopic Survey (BOSS, Dawson et al. 2013).BOSS was a survey conducted as part of the third stage of the larger Sloan Digital Sky Survey (SDSS, York et al. 2000), which collected optical spectra from more than 1.5 million targets using the 2.5-m Sloan Telescope (Gunn et al. 2006) at Apache Point, New Mexico.BOSS covered roughly 10,000 deg 2 of the sky in two hemispheres, referred to as the North and the South Galactic caps (NGC and SGC, respectively).
Our analysis is focused on the CMASS galaxy sample, which is dominated by luminous red galaxies (LRG) that were selected on SDSS multicolour photometric observations (Gunn et al. 1998;Gunn et al. 2006).CMASS is nearly complete down to stellar mass of  * ≈ 10 11.3 M ⊙ for  > 0.45 (Maraston et al. 2013), and covers a redshift range 0.4 ≲  ≲ 0.7.For this paper, we impose a more stringent redshift cut 0.45 ≤  ≤ 0.6 to avoid regions where the galaxy number density drops abruptly, which can potentially bias our model predictions.Additionally, we restrict the analysis to the NGC for simplicity.After imposing these restrictions, this results in a sample with a total volume of ≈ 1.4 ( ℎ −1 Gpc) 3 , an effective volume of ≈ 1.1 ( ℎ −1 Gpc) 3 , and an average number density of ≈ 3.5 × 10 −4 ( ℎMpc −1 ) 3 .
We use the DR12 large-scale structure catalogues provided by the BOSS collaboration 1 (Reid et al. 2016).These catalogues include angles and redshifts for each galaxy, which we convert to comoving Cartesian coordinates by adopting a flat-ΛCDM fiducial cosmology characterized by a matter density parameter Ω m = 0.315, which closely matches the Planck 2018 best-fit cosmology, assuming base-ΛCDM (Planck Collaboration et al. 2020).The BOSS collaboration also provides a set of random catalogues that follow the footprint and radial selection of CMASS galaxies, but with no intrinsic clustering, which are used to estimate the overdensity field as described in the following sections.

Two-point clustering
Galaxy clustering is usually characterized in terms of the two-point correlation function (2PCF),  gg (), which, in its simplest form, quantifies the excess probability d of finding a galaxy in a volume d, separated at a distance  from another galaxy, with respect to an unclustered Poisson distribution: where   is the mean number density of galaxies in the sample.In the presence of redshift-space distortions (RSD, Jackson 1972;Kaiser 1987) or Alcock-Paczynski distortions (AP, Alcock & Paczynski 1979), the galaxy distribution appears anisotropic to the observer.To capture this anisotropy, a convenient choice is to bin the correlation function in terms of  and , where  is the redshift-space pair separation and  is the cosine of the angle between the vector connecting the two galaxies and the observer's line of sight.
A number of estimators have been proposed to measure the 2PCF from observational data.A robust and commonly used estimator is the one introduced by Landy & Szalay (1993), where GG is the normalized number of galaxy pairs in the (, ) bin, while GR and RR are the normalized galaxy-random and randomrandom pairs, which make use of the unclustered random catalogues.
It is useful to separate out the different angular components of the (, ) correlation function by decomposing it into multipole moments, defined by with P ℓ the Legendre polynomials.We measure  gg (, ) in CMASS using pycorr2 , which is a wrapper around a modified version of the CorrFunc pair-counting code (Sinha & Garrison 2020).We focus the analysis on the monopole ( 0 ) and quadrupole ( 2 ) moments.The correlation functions are measured in 241  bins from −1 to 1, and radial bins with scale-dependent widths: 1ℎ −1 Mpc bins for  ∈ [0, 4] ℎ −1 Mpc, 3ℎ −1 Mpc bins for  ∈ (4, 30] ℎ −1 Mpc, and 5ℎ −1 Mpc bins for  ∈ (30, 151] ℎ −1 Mpc.Galaxies are appropriately weighted during the pair counting to consider various observational systematics that can bias the clustering measurements.The total systematic weight for each galaxy is given by where  sys ,  fc , and  zf account for imaging systematics, fiber collisions, and redshift failures (Ross et al. 2016).This is additionally multiplied by a weight  FKP , which optimally weights the contribution of galaxies based on their redshift-dependent number density (Feldman et al. 1994), where  0 = 10 4 ℎ3 Mpc −3 .The total weights for each galaxy are then given by while points from the random catalogue are only weighted by  FKP .

Density-split clustering
The density-split clustering (DSC) method (Paillas et al. 2021) characterizes galaxy clustering in bins of environmental density, with the aim of extracting and combining the cosmological information coming from distinct environments of the cosmic web.We apply the density-split algorithm to the CMASS galaxy sample using our publicly available code 3 , with slight modifications to the algorithm presented in Paillas et al. (2023) to account for the non-uniform survey geometry.We start by painting the CMASS galaxies and randoms to a rectangular grid that fully encompasses the survey volume, and we estimate the overdensity field as where G and R are the normalized galaxy and random counts in each cell, weighted as given in Eq. ( 4).We smooth the overdensity field with a Gaussian filter of radius   = 10 ℎ −1 Mpc, and then we sample it using cloud-in-cell interpolation at  query query positions, which are taken from the CMASS random catalogue.Here we set  query to 5 times the number of galaxies in the catalogue, and split those query points into 5 quintiles, according to the overdensity at each location.Figure . 1 shows the probability distribution function (PDF) of galaxy overdensity measured at the query positions.The overdensity field shows a non-Gaussian PDF with significant skewness and kurtosis.This shape is a consequence of the growth of structure being bounded at Δ = −1 from below (regions completely devoid of galaxies), while no such constraint is present at the positive Δ end.The distribution peaks at negative overdensities, reflecting that the average region in the Universe is underdense due to the larger volume occupied by voids.The division into quintiles is demarcated by the different colours in the figure.We label the quintiles as Q  , where  goes from 0 to 4 from lower to higher densities.In what follows, we discard Q 2 from the clustering analysis, since all five quintiles are not independent from each other. 4igure 2 shows the comoving number density of quintile positions as a function of redshift, which is observed to closely follow that of the galaxies.This is given by construction, as the query positions are sampled from the clustering random catalogues, which in turn are constructed to match the footprint and radial selection of the galaxies.
Once the density quintiles are defined, we estimate the quintilegalaxy cross-correlation function using the estimator (Landy & Szalay 1993) where QG, QR, RG, RR are the normalized quintile-galaxy, quintilerandom, galaxy-random, and random-random pair counts.We note that this assumes that the random catalogue is the same for galaxies and quintiles, which is justified by the good match between the () distributions from Fig. 2. The autocorrelation functions of DS quintiles are estimated as with QQ being the normalized quintile-quintile pair counts.We adopt the same binning scheme as for the galaxy 2PCF, as described in Sect.2.2.1, and we weight the galaxy pairs according to Eq. (6).
Figure 3 shows the measured multipoles from DSC and the galaxy 2PCF.The DSC auto-and cross-correlation functions cover a broad range of amplitudes, as the quantiles trace the underlying matter density field in different ways.Q 0 & Q 1 are underdense regions with negative linear bias parameters that usually range from -3 to -1, whereas Q 3 & Q 4 tracer overdensities and have positive bias parameters that can range from 1 to 3 (Paillas et al. 2021).Two signature features are spotted in all profiles: a transition regime around 25 ℎ −1 Mpc, which is due to the scale that was used to smooth the overdensity field and define the quintiles, and the peak/valley around 100 ℎ −1 Mpc, which is an imprint of the BAO that originated from sound waves in the photon-baryon plasma prior to recombination.
The quadrupoles show a large degree of anisotropy in DSC.Two main factors contribute to this anisotropy.Firstly, there is an RSD effect caused by the dynamics of galaxies around different density environments.Kaiser-like motions on large scales and random motions on small scales (Fingers of God) cause distinct RSD patterns on the clustering of each quintile, similar to the well-known effects seen in the galaxy 2PCF.Secondly, there is an RSD effect imprinted on the quintile positions themselves, which is the product of identifying the density quintiles in redshift space.Paillas et al. (2023) showed that splitting densities in redshift space causes selection effects that induce distortions in the clustering of the quitiles, which manifests itself as a quadrupole moment in the DSC autocorrelation functions.This is similar in nature to the selection effect that is produced when cosmic voids are identified in redshift space (Chuang et al. 2017;Nadathur et al. 2019;Correa et al. 2020).
We reserve the discussion of the model fits (solid lines) for Sect.4.1.

Mock galaxy catalogues
Throughout this work, we use different types of mock galaxy catalogues to study various aspects of the analysis.The MD-Patchy mocks are used to estimate the sample variance associated to our observed clustering measurements.The Nseries mocks are used to validate the theory model we apply to CMASS.The AbacusSummit mocks are used to train our simulation-based model for galaxy clustering, as presented in our companion paper (Cuesta-Lazaro et al. 2023).Each of these mocks are described in more detail below.

MD-Patchy mocks
To estimate the covariance of the data vector, we use the MultiDark-Patchy mocks (MD-Patchy, Kitaura et al. 2016), a suite of 2048 mock galaxy catalogues that were designed to match the footprint, redshift distribution and halo occupation distribution of the BOSS DR12 galaxy samples.The mocks are based on approximate lightcones generated with augmented Lagrangian perturbation theory (ALPT) using the patchy code (Kitaura et al. 2014(Kitaura et al. , 2015)).ALPT is based on a combination of second-order Lagrangian perturbation theory on large scales and the spherical collapse model on smaller scales.patchy populates dark matter haloes using a subhalo abundance matching prescription that is calibrated from N-body simulations from the Big MultiDark Suite (Klypin et al. 2016)  8 = 0.8288, a tilt of the primordial power spectrum   = 0.9611, and a dimensionless Hubble parameter ℎ = 0.6777.

Nseries mocks
To study systematic errors in our theory model, we use the Nseries cutsky mocks, a collection of 84 mock catalogues that were designed to match the clustering, footprint and radial selection of the CMASS galaxy sample.The cutsky mocks are constructed from a base set of seven independent full N-body simulations, with a box side length of 2.6 ℎ −1 Gpc and a mass resolution of 1.5 × 10 11 ℎ −1 M ⊙ .Dark matter haloes are populated with a halo occupation distribution prescription, with parameters chosen to best model the clustering of the BOSS DR12 CMASS sample.Each box is trimmed and rotated in different ways to produce 12 cutsky mocks that match the geometry of the CMASS sample, which results in a total of 84 pseudo-independent cutsky catalogues.The mocks are then passed through the same fiber assignment code that was used for BOSS, so that they faithfully reproduce the angular variations of fiber collisions in the data (Hahn et al. 2017).

AbacusSummit mocks
To train and validate our clustering emulators, we use Abacus-Summit, a suite of cosmological N-body simulations (Maksimova et al. 2021) designed to meet the simulation requirements of the Dark Energy Spectroscopic Instrument (Levi et al. 2019).The simulations were run with the Abacus N-body code (Garrison et al. 2019(Garrison et al. , 2021)), comprising different volumes, resolutions, and cosmologies.The base simulations follow the evolution of 6912 3 dark matter particles in a (2 ℎ −1 Gpc) 3 volume, with a mass resolution of 2×10 9 ℎ −1 M ⊙ .There are 97 cosmology variations in total, exploring an eight-dimensional parameter space around the Planck18 primary ΛCDM cosmology (PL18, Planck Collaboration et al. 2020): Here,  cdm = Ω c ℎ 2 and  b = Ω b ℎ 2 are the physical cold dark matter and baryon densities, d  /d ln  is the running of the spectral tilt,  eff is the effective number of ultra-relativistic species,  0 is the present-day dark energy equation of state, and   captures the time evolution of the dark energy equation of state.The simulations assume a flat spatial curvature, and the dimensionless Hubble parameter ℎ is calibrated to match the Cosmic Microwave Background (CMB) acoustic scale  * to the PL18 measurement.
For the training and validation of the emulators, we restrict to the Table 1.List of cosmological and HOD parameters used in our analysis.For each, we quote the parameter symbol, the prior distribution, the fixed value in the baseline model (where appropriate), and the physical interpretation.We note that the prior distribution for all parameters is uniform, with the exception of the baryon density, for which we adopt a normal distribution with a mean and dispersion as specified.c001-004 Secondary cosmologies, including a low  cdm choice (WMAP7, Komatsu et al. 2011), a CDM choice, a high  eff choice, and a low  8 choice.
c100-126 Linear derivative grid providing paired simulations with small negative and positive steps in the eight-dimensional cosmological parameter space from Eq. ( 10).
c130-181 An emulator grid around the c000 cosmology that provides a wider coverage of the cosmological parameter space.
In addition to the base simulations, there are multiple realizations of smaller boxes with a side length of 500 ℎ −1 Mpc at the c000 cosmology, which can be used for covariance estimation.Throughout the rest of the paper, we will refer to these simulations as AbacusSmall, and to the base simulations simply as AbacusSummit.
Group finding is done on the fly, using a hybrid Friendsof-Friends/Spherical Overdensity algorithm, dubbed CompaSO (Hadzhiyska et al. 2021).As described in Sect.3.3, we populate these halo catalogues with galaxies using a halo occupation distribution prescription, which are then used to obtain the clustering measurements for our training data.

Galaxy-halo connection model
In the current paradigm of cosmology, galaxies are thought to form and evolve within dark matter halos, which are large structures that form as a result of the gravitational collapse of overdensities in the Universe.The halo occupation distribution (HOD) is a statistical model that describes how galaxies are distributed within dark matter halos.
A well-suited model for LRG is the base halo model from Zheng et al. (2007), where the average number of central galaxies in a halo of mass  is given by where erf () denotes the error function,  cut is the minimum mass required to host a central, and  is the slope of the transition between having zero to one central galaxies.The average number of satellite galaxies is in turn given by where  cut gives the minimum mass required to host a satellite,  1 is the typical mass that hosts one satellite, and  is the power law index for the number of galaxies.
We use the AbacusHOD package, which is highly efficient and contains a wide range of HOD extensions (Yuan et al. 2022a).For this work, we extend the base model to modulate galaxy peculiar velocities via the parameters  vel,c , which parametrizes the velocity bias between the central galaxy and the halo centre, and  vel,s , which parametrizes the velocity bias between the satellite galaxies and the local dark matter particles.When no velocity bias is present,  vel,c = 0 and  vel,s = 1, in which case centrals perfectly follow the velocity of halo centres, and satellites perfectly match the velocity of the underlying dark matter particles.
We introduce two additional parameters to account for galaxy assembly bias:  cen and  sat , which add environment-based secondary bias for centrals and satellites, respectively.Here, the environment is defined as the smoothed matter density around the halo centres, using a top-hat filter of radius   = 5 ℎ −1 Mpc.When no secondary bias is present,  cen =  sat = 0. Positive/negative values of these parameters indicate a preference for galaxies to form in haloes around less/more dense environments, respectively.

Clustering emulators
In the context of galaxy clustering and cosmology, an emulator refers to a computational model or algorithm that approximates the predictions of expensive or time-consuming simulations or calculations.They are used to efficiently and accurately generate predictions from cosmological models without resorting to additional simulations.
Emulators are trained on a set of pre-computed simulations, often referred to as a training or calibration set, which cover a wide range of model parameter values and capture the desired properties of interest of an observable.As such, the emulator can learn how this observable responds to changes in cosmology or galaxy-halo connection parameters.Once the emulator is trained, it can rapidly predict the desired model predictions for any given set of input parameters.
In Cuesta-Lazaro et al. ( 2023), we present our galaxy 2PCF and DSC emulators, which are based on HOD catalogues constructed from the AbacusSummit simulations.Here, we present a brief description of how the emulators are constructed.
We start from the dark matter halo catalogues from the Abacus-Summit snapshots at  = 0.5, spanning 85 different cosmologies within the eight-dimensional  0   ΛCDM parameter space defined in Eq. ( 10).Using the abacushod code (Yuan et al. 2022a), we populate dark matter haloes with a nine-parameter extended HOD framework (Sect.3.2), generating 100 unique HOD variations per cosmology, where the HOD parameters are sampled from a Latin Hypercube to ensure an optimal sampling of the parameter space.When the number density of an HOD catalogue is above the average number density of the CMASS sample ( gal ≈ 3.5 × 10 −4 (ℎ/Mpc) 3 ), we invoke an incompleteness parameter  ic to downsample the catalogue down to the target number density.Under the distant-observer approximation, we map the positions of galaxies to redshift space by perturbing their positions with their peculiar velocities along one cartesian axis of the simulation box chosen to represent the line of sight.We repeat this procedure for each of the three axes of the simulation boxes, effectively generating three pseudo-independent redshift-space catalogues for each HOD variation, from which we can average the clustering measurements to reduce cosmic variance later on.
For each of these mock catalogues, we run the DSC pipeline, using a mesh resolution of  cell = 5 ℎ −1 Mpc, a number of random query points  query equal to five times the number of galaxies in each catalogue, a smoothing radius for the Gaussian filter of   = 10 ℎ −1 Mpc, and five density quintiles.We compute the galaxy 2PCF and the density-split correlation functions in bins of  and , and we decompose them into their multipole moments.
We split the HOD catalogues into training, validation and test sets, and for each clustering statistic we construct separate fully-connected neural networks, which take the cosmological and HOD parameters as an input, and return the monopole and quadrupole moments of the correlation functions.For training and validation, we use cosmologies c100-126 and c130-181 , whereas the rest are reserved for the test set.The hyperparameters of the neural networks are calibrated to minimize the validation loss.Overall, we observe that the emulators produce model predictions with percent-level accuracy for the full range of scales, when tested against the test simulation boxes.

Likelihood
When fitting the emulator to the CMASS data, we define the loglikelihood as where  data is the observed data vector,  model is the emulator prediction, and C is the covariance matrix, which includes three contributions to the error budget: Here, C data is the term associated with the sample variance of the data vector, which is estimated from multiple realizations of the MD-Patchy mocks: where  patchy = 2048.The simulations used for training are at a fixed phase, which could be different from the true underlying phase of the Universe.In other words, the cosmic variance is frozen in our emulator predictions.To account for this, we add an extra sample variance contribution to the error budget, associated with the finite size of the training simulations.We estimate this covariance using multiple mock realizations of the AbacusSummit fiducial cosmology with a fixed set of HOD parameters with high likelihood: where  sim = 1800.C emu accounts for the intrinsic error in the model predictions due to an imperfect emulation.This term is calculated by computing a covariance matrix from the difference between the emulator predictions and measurements from a set of test simulations with known cosmologies and HOD parameters, Δ: where the overline denotes the mean across all test simulations, and  test = 600.When calculating the likelihood, both C data and C sim are multiplied by a factor P before inversion (Percival et al. 2022): where Here,   is the number of simulations used to estimate the covariance,   is the length of the data vector and   is the number of parameters that are being fitted.This corrects for the fact that even though Eq. ( 16) & ( 17) are unbiased estimates of the covariance matrices, the inversion leads to biased parameter constraints.We sample the posterior distribution of parameters using the dynesty nested sampler (Speagle 2020), which also provides an estimation of the Bayesian evidence.As listed in Table 1, we assume flat priors for all parameters except for the baryon density, for which we adopt a BBN-informed Gaussian prior (Aver et al. 2015;Cooke et al. 2018;Schöneberg et al. 2019) for our baseline analysis: In later sections, we also explore using a flat prior range for  b , finding that it has little impact on our conclusions and the reported constraints for other parameters.

𝜃 * prior
In the simulations used to train our emulator, the value of the dimensionless Hubble parameter ℎ was chosen to fix the CMB acoustic scale  * (100 * = 1.041533) matching the best-fit measurement from PL18 (Planck Collaboration et al. 2020;Maksimova et al. 2021).This means that effectively we apply a fixed  * constraint as a prior on our emulator results.This ensures that we only consider models within the part of parameter space where the emulator has been trained.As a result, ℎ is not a free parameter in our model, but is determined by the sampled parameters and the  * constraint.We use class (Blas et al. 2011) to derive ℎ at each point in our chains, and we then use this to also obtain Ω m , calculated as where   = 0.00064420 accounts for 60 meV neutrinos (also the choice of the base PL18 cosmology), which is always fixed in our model.
A consequence of this prior on  * is that the parameter constraints we obtain below do not come exclusively from late-Universe clustering measurements but rather from a combination of galaxy clustering and information on the CMB acoustic scale, although they do not use

RESULTS
In this section, we apply our emulator framework to the CMASS clustering data and present our main cosmological constraints, along with tests for model systematics.We focus on the constraints for cosmological parameters and reserve the discussion about HOD constraints for Appendix A.

Model fits
The solid lines in Fig. 3 show the best-fit base-ΛCDM model to the measured multipoles, with the  2 , likelihood, and Bayesian evidence values reported in Table 2.The  2 per degree of freedom for the galaxy 2PCF + DSC combination is 1.66, while for the 2PCF-only fit it is 1.06, which is comparable to the reported values for the smallscale 2PCF fit from Yuan et al. (2022b) and the configuration-space BAO fit from Ross et al. (2016).
We observe some hints of the models under-predicting (or overpredicting, in the case of negative density contrasts) the observed clustering at scales  > 90 ℎ −1 Mpc.Such excess large-scale clustering has been observed before in the galaxy 2PCF analysis of BOSS data.In Ross et al. (2016), the authors find that the observed monopole shows an apparent excess with respect to the mean of the MD-Patchy mock catalogues, although it is argued that the mismatch is of low statistical significance.In Satpathy et al. (2017), the authors fit the correlation function using a model based on Convolutional Lagrangian Perturbation Theory, also finding that the model underpredicts the monopole data at large scales.Here, we observe a similar trend when comparing the measured multipoles to our emulator predictions, although the level of discrepancy is never greater than two standard deviations for all scales considered, irrespective of the summary statistic.As informed by the  2 values, the quality of the fit is better than might be guessed by the eye, since the separation bins are correlated.In Sect.4.4, we show that the mean clustering signal of the Nseries mocks monopole is also lower than the CMASS data, although still consistent within 2.We also observe that the emulator is able to fit the mean of the mocks with much higher accuracy than the data, which we deem reasonable given that our model was trained to make predictions for the ensemble average of the clustering statistics.
Although some of the observed offsets between our best-fit model and CMASS could be partially attributed to statistical fluctuations, we should bear in mind the possibility that there are residual systematic effects in the clustering data that are not fully corrected by the weighting procedure.Lavaux et al. (2019) performed a field-level inference of BOSS data, finding evidence of residual systematics on the spectroscopic data that so far remain unexplained, which can induce correlated modulations of the order of 30 per cent on the sky for CMASS.In addition, it is worth noting that the systematic weights for the large-scale structure catalogues provided by the BOSS collaboration were mainly validated for two-point statistics, and it is possible that they are sub-optimal for alternative clustering methods such as DSC.We plan to explore this topic further in future work.

Base-ΛCDM constraints
In this section, we explore the constraining power of DSC and the galaxy 2PCF on the base ΛCDM parameters, letting  b ,  cdm ,  8 ,   , and the HOD parameters vary during the fit.To generate the model predictions, we query the emulator fixing the remaining cosmological parameters to their baseline values  0 = −1,   = 0, d  /d ln  = 0, and  eff = 3.0146.
Figure 4 shows the 2D posterior distributions on ΛCDM, marginalized over the HOD parameters.We also over-plot the base-ΛCDM constraints from the Planck_TTTEEE_lowl_lowE likelihood from PL18 to facilitate comparison.The reported best-fit values, along with the means of the marginalized posteriors and their dispersion are listed in Table 3.It is worth noting that the best-fit values can be shifted with respect of the mean of the posterior due to its non-Gaussian shape.
The galaxy 2PCF by itself constrains  cdm ,  8 , and   with a 2.7, 6.2, and 3.2 per cent precision, respectively.Adding the DSC multipoles tightens the constraining power to a precision of 1.7, 3.8 and 1.8, respectively.The constraints for baryon density are entirely dominated by the BBN prior.Overall, both posterior distributions are largely consistent with the PL18 results, with a 0.04 and 0.6 difference in the mean values of  cdm and  8 between our baseline fit and PL18.
As described in Sect.3.5, we use the  * constraint to obtain ℎ and Ω m as derived parameters from our chains, and also show the results for them in Fig. 4. We obtain ℎ = 0.6793 ± 0.0070 and Ω m = 0.3122 +0.0094 −0.011 when fixing to the base ΛCDM model.As Fig. 4 shows, we do not observe any significant degeneracy between Ω m and the baryon density (and our constraints on  b are also prior-dominated).
This means that, to a good approximation,  * ∝ Ω m ℎ 3.4 , the expected behaviour for fixed Ω b ℎ 2 within flat-ΛCDM models (Percival et al. 2002), and thus DSC and the galaxy 2PCF results both describe this narrow degeneracy direction in the Ω m -ℎ plane.This is close to but not the same as the similarly tight degeneracy obtained from the the CMB by PL18, which corresponds to constant Ω m ℎ 3 .The degeneracy direction for models with the same  * depends weakly on Ω b ℎ 2 through the sound horizon.PL18 uses the CMB data to also constrain Ω b ℎ 2 , leading to an anti-correlation between constraints on Ω b and Ω m , which then results in a different degeneracy in the Ω m -ℎ plane.
Within ΛCDM, the linear growth rate of structure can be approximated as We use this expression to derive a value for   8 at the effective redshift of CMASS from our chains 5 .We obtain   8 ( = 0.525) = 0.462±0.020,which is a 4.3 per cent constraint from the combination of DSC and the galaxy 2PCF.We contrast this constraint with the PL18 prediction and other clustering studies in Fig. 5. Our constraint is slightly lower than the mean of PL18, but consistent at the 1 level.We also find good agreement with the results of the main BOSS clustering analysis (Alam et al. 2017), but we get a factor of 1.9 better precision using the DSC + galaxy 2PCF combination, compared to their middle redshift bin that covers 0.4 <  < 0.6.We find a higher mean   8 than Yuan et al. (2022b) (from hereon Y22), although there is still consistency at the 0.7 level.Interestingly, Y22 constrain   8 to a precision similar to our baseline analysis, even though they only use the galaxy 2PCF at scales < 30 ℎ −1 Mpc.Our 2PCF-only fit is a factor of 1.8 worse than theirs, which seems puzzling given that their model is based on the same set of simulations used in this study.The main reason we attribute this difference to is that while we estimate the intrinsic emulator error using a test set of simulations that covers the whole prior range in the cosmology & HOD parameter space, Y22 estimate the error using only HOD catalogues that have a high likelihood with respect to the measured data vector.Our emulator error will generally tend to be more conservative, as it also considers the error around regions near the edge of the priors, which are usually harder to emulate.It also dominates the error budget at small scales, precisely where Y22 extract most of their information.Effectively, this means we are discarding some small-scale information which they keep.We have verified that by artificially lowering or removing the emulator error, the precision of our galaxy 2PCF constraints closely matches that of Y22.Other reasons that could potentially add to this are that Y22 use Gaussian priors on the HOD parameters while we use flat priors, and that Y22 estimate the data covariance through jackknife resampling, while our covariance is estimated from mock catalogues.
Our results are 2.2 higher than those reported by Zhai et al. (2023), who find   8 ( = 0.55) = 0.396 ± 0.022.Zhai et al. (2023) train an emulator using the Aemulus suite of simulations to fit the small-scale clustering of BOSS galaxies using the galaxy 2PCF between 0.1 ℎ −1 Mpc and 60 ℎ −1 Mpc.Although we have not tested our emulator against Aemulus directly, our companion paper (Cuesta-Lazaro et al. 2023) shows that our model can recover unbiased cosmological constraints when fitted to mock galaxy catalogues generated with a different N-body code and galaxy-halo connection model, 5 As we are doing a direct fit of ΛCDM parameters and marginalizing over ℎ, our results are not affected by the problems of the standard template-based clustering analyses pointed out in Sánchez (2020).2023) use Eulerian perturbation theory in combination with a halo model calibrated on N-body simulations to model the power spectrum multipoles from BOSS DR12 up to  = 0.2 ℎMpc −1 , finding   8 ( = 0.61) = 0.455 ± 0.026, which is in excellent agreement with our constraint, albeit being at an effective redshift slightly higher than our sample.
Comparing with a full-shape analysis based on the Effective Field Theory of Large-Scale Structures (EFTofLSS, Carrasco et al. 2012), we find an   8 value than is 1.7 higher than the one derived by d'Amico et al. ( 2020), who fit the power spectrum multipoles of BOSS DR12 galaxies up to  = 0.2 ℎMpc −1 and report   8 ( = 0.55) = 0.399 ± 0.031, which is closer to the Zhai et al. (2023) estimate.Although our Ω m constraints agree to within 0.1 with d'Amico et al. ( 2020), their predicted amplitude of the primordial power spectrum,   , which is 2.3 lower than PL18, drives their derived   8 to lower mean values.
Finally, Semenaite et al. (2022) presented a clustering analysis of BOSS and eBOSS data, using information of the full shape of the BOSS clustering wedges presented by Sánchez et al. (2017), and the multipoles from eBOSS quasars from Hou et al. (2021).They derive  8 = 0.815 ± 0.044 and Ω m = 0.290 +0.012 −0.014 , which agree with our constraints at the 0.4 and 1.2 level, respectively.

Extended-ΛCDM constraints
To explore potential deviations from ΛCDM, we have analyzed a grid of well-motivated extensions to the base model.
Figure 6 and Table 4 summarize the results of single-parameter extensions to the base ΛCDM model.All these fits have been run using the baseline configuration of data vectors, scale cuts and model prescription, except for the addition of a single parameter to the cosmological model, which are incorporated one at a time.We do not find compelling proof supporting any of these extensions, as the marginalized posteriors for the additional parameters generally overlap with the base model within one standard deviation.

Running of the spectral index
The first extension we look at is in regard to the scale dependence of the primordial density fluctuations.Here we characterize the primordial power spectrum as a power law with a normalization amplitude   , a spectral index   and its first derivative with respect to ln   (Yuan et al. 2022b;Kobayashi et al. 2022;Zhai et al. 2023), an EFTofLSS analysis (d'Amico et al. 2020), and a halo perturbation theory analysis (Yu et al. 2023).In some cases, the redshifts of the measurements have been slightly shifted horizontally for visual clarity.
(also known as the running): where  0 is a pivot wavenumber used to specify the point at which the power spectrum is normalized.Our base-ΛCDM constraints from Sect. 3 have assumed a zero running for the spectral index, finding   = 0.970 ± 0.018.CMB experiments also favour   < 1, which is predicted by common single-field slow-roll inflationary models (Mukhanov 2007).Such models also predict a very small running, since it is only second-order in inflationary slow-roll parameters (Kosowsky & Turner 1995), but it is possible to construct valid models that predict larger values.We find d  /d ln  = −0.005± 0.013, which is consistent with zero running.The 1D marginalized posterior is shown in the left-hand side panel of Fig. 6, where we also overplot the constrain from the Planck_TTTEEE_lowl_lowE likelihood from PL18, who find −0.0055± 0.0067.We warn the reader that the parameter range for the running shown in Fig. 6 corresponds to our full prior range, beyond which we do not have simulations to sample the parameter space with our model.The posterior distribution quickly approaches zero near the prior walls, giving us confidence that our 1 limits are not prior-dominated.However, we cannot confidently provide 3 limits given this restriction.

Effective number of relativistic species
A relic neutrino background is a generic prediction of the standard hot Big Bang model (Lesgourgues & Pastor 2014).The constraints on the properties of relic neutrinos and other relativistic species beyond the Standard Model of particle physics is of special interest for large scale structure analyses.The combination of CMB experiments, galaxy and supernova surveys have put the tightest upper limits on the sum of neutrino masses (Emas & Wulandari 2019), whereas Planck has constrained the density of light relics with sub-percent level precision (Planck Collaboration et al. 2020).In the instantaneous neutrino decoupling limit, the density of radiation in the Universe (besides photons) can be written as (Lesgourgues & Pastor 2014): where  eff = 3, usually called the effective number of relativistic species, is a convenient parametrization of the relativistic energy density of the Universe beyond just photons, in units of a single neutrino.Detailed calculations that go beyond the instantaneous neutrino decoupling limit, including neutrino oscillations, predict  eff ≈ 3.046 (de Salas & Pastor 2016).
For the base-ΛCDM constraints from Sect. 3, we fixed  eff to the baseline value of 3.046.In this section, we vary this parameter during the inference analysis, finding  eff = 3.02 +0.24 −0.27 , which is a 8.4 per cent precision constraint, in excellent agreement with the PL18 measurement.

Dark energy equation of state
One of the main goals of modern observational cosmology is elucidating the nature of the accelerated expansion of the Universe.In the base-ΛCDM model, dark energy takes the form of a cosmological constant that has an equation of state  0 ≡ / = −1, where  and  represent the pressure and density of the fluid.In this section, we explore models with a constant  by letting it vary as a free parameter in the fit.
The right-hand side panel of Fig. 6 shows the marginalized posterior of  for the combination of the galaxy 2PCF and DSC.We find  0 = −0.956constant at the 1 level.CMB data alone (blue contour) does not put a very tight constraint on  0 , as it is a  ≈ 1100 measurement.The grey contours show the posterior resulting from the combination of the Planck_TTTEEE_lowl_lowE PL18 likelihood with late-time probes of the expansion rate, including BAO measurements from BOSS DR12 (Alam et al. 2017), SDSS-MGS (Ross et al. 2015), and 6dFGS (Beutler et al. 2011), as well as type Ia supernovae distance measurements from Pantheon sample (Scolnic et al. 2015), and local estimates of the Hubble parameter from Milky Way Cepheid variables from Riess et al. (2018).This combination tightens the marginalized posterior to  0 = −1.041+0.060 −0.053 , which agrees with our constraint at the 1.1 level.
We emphasize that for our results from galaxy clustering, we have adopted the prior constraint for the acoustic scale  * from PL18, so our constraints on  0 are not fully independent from the Planck_TTTEEE_lowl_lowE likelihood.

Test for systematics
In this section, we carry out tests on mock galaxy catalalogues to look for systematics in the theorerical modelling, and we assess the robustness of our cosmological constraints to different choices in our inference pipeline.

Recovery tests on the Nseries mocks
We begin by testing our pipeline on the Nseries mock galaxy catalogues, which were calibrated onto the clustering of the CMASS NGC galaxy sample, matching its footprint and radial selection (Sect.3.1.2).We measure the galaxy 2PCF and density-split multipoles from each of the 84 mock realizations, and analyze each mock using the baseline configuration of our pipeline as we did in Sect.3. We estimate the covariance matrix of the data vectors from the MD-Patchy mocks, following the same procedure as described in Eq. ( 16).
The left panel of Fig. 8 shows the distribution of recovered bestfit values from the 84 fits.The true cosmology of the simulations, which is shown by the vertical red lines, is well within the 68 per cent confidence region of the distribution, showing that our clustering pipeline is able to recover unbiased cosmological constraints even in the presence of complex survey masks, fiber collisions and nonuniform radial selection functions.
As a complementary test, we proceed to average the data vectors across the 84 mock realizations, and perform cosmological fits on the mean data vectors.The right panel of Fig. 8 shows the recovered cosmological parameters in this setup, using two different covariance matrices.The grey contours show results with the usual covariance matched to the volume covered by the CMASS sample used throughout the paper.The black contours show results where the covariance is rescaled to match the comoving volume covered by 84 realizations6 of the Nseries suite, which amounts to roughly 120 (ℎ −3 Gpc 3 ), after applying the redshift cuts to match our data sample, 0.45 <  < 0.6.Even using this large volume, the true cosmology of the simulations falls within one standard deviation of the marginalized parameter posterior distributions, highlighting the robustness of our pipeline even for datasets far larger than the one analysed in this paper.
Figure 10 shows the best fit to the galaxy 2PCF measured from the mean of the Nseries samples (using the covariance associated to a single CMASS volume), where we also overplot the measurement from the CMASS sample for comparison.Overall, the fit to the Nseries mocks is accurate at all scales, with the best-fit model always falling within one standard deviation of the error bars.This can be contrasted with the fit to the CMASS data seen earlier in Sect.4.1, where the data measurement shows a lower and higher clustering amplitude at intermediate and large scales compared to the model, respectively.CMASS also shows a similar difference in clustering with respect to the mean of the Nseries mocks.The  2 between the mocks and data monopole is 30 for 23 degrees of freedom, which corresponds roughly to a 2 shift.Although we do not show the other statistics for brevity, we have observed the same trends for DSC.However, it is worth keeping in mind that the cosmology of the Nseries simulations differs substantially from the fiducial cosmology we adopted to convert redshifts to distances, so the Nseries multipoles could be affected by more severe AP distortions than the CMASS data, if the true cosmology of the Universe is closer to our fiducial cosmology compared to the mocks.
Figure 7 shows the distribution of log-evidence from the fits to the 84 Nseries mocks and CMASS.We include the results obtained using two different scale cuts:  min = 1 ℎ −1 Mpc and  min = 50 ℎ −1 Mpc.We see that regardless of the scale cut that is used, the evidence of the CMASS fit is significantly lower than for any of the Nseries mocks.This means that, overall, the model is a much better description of the mocks than it is of the data, even when marginalizing over all log-evidence .Logarithm of the Bayesian evidence obtained from individual fits to the Nseries mocks (histograms) and to the BOSS DR12 CMASS data (dashed lines).We show results for two different minimum scale cuts:  min = 1 ℎ −1 Mpc (green) and  min = 50 ℎ −1 Mpc (pink).The evidence of the CMASS fit is significantly lower than any of the Nseries mocks in both cases, but this discrepancy is reduced when excluding small scales.
the parameter space of the model.This supports the hypothesis that there might be residual systematic effects in the CMASS clustering catalogues that are currently unaccounted for by the weighting procedure (Lavaux et al. 2019).While removing all the information below  = 50 ℎ −1 Mpc reduces the gap between Nseries and CMASS, it is still statistically unlikely that the large offsets between the baseline model and the data can be fully explained by a random fluctuation due to sample variance from the observations, if we assume that the mocks are a faithful reproduction of the CMASS data.This could be a combination of the above-mentioned systematic effects with insufficient models of the halo-galaxy connection on small scales.
Given the small number of Nseries mocks, we cannot translate these findings into more quantitative statements.

Scale cuts
Another important aspect we study with the Nseries mocks is the choice of scales that are used in the main data analysis.Figure 9 shows the marginalized constraints on  cdm ,  8 ,   , and   8 , adopting eight different minimum scale cuts, ranging from 1 ℎ −1 Mpc to 60 ℎ −1 Mpc.For these measurements, we use the mean of the 84 Nseries mocks as the data vector, and a covariance associated to a single CMASS volume.We observe that for Nseries, the mean values of the marginalized posteriors are very stable when changing the scale cuts for all parameters that are considered.It is interesting to note that the size of the error bars does not shrink beyond ≈ 20 ℎ −1 Mpc.There are two main reasons for this.The first one is the inclusion of the model uncertainty in the covariance matrix that is used for the calculation of the likelihood [Eq.( 15)].At these scales, although the emulator error has a percent-level accuracy, it starts to dominate the total error budget for the monopole of the 2PCF and the density-split CCF.This ensures that the cosmological inference is always robust, even when including scales where the emulator is less accurate than the precision of the data, which comes at the expense of not being able to extract all the information there is available.The second reason is that even though we are imposing a minimum scale cut in the multipoles, the DSC quintiles are defined using small-scale informa-tion from the density field, also propagates into the multipoles at larger scales, as shown in Paillas et al. (2023).Figure 9 also shows how the constraints from CMASS data change depending on the minimum scale cut.Although the results we get using  min = 1 ℎ −1 Mpc are consistent with those obtained with more conservative cuts to within 1, we observe some interesting trends with scale. cdm shows some tendency towards larger mean values the more small-scale information we include.Furthermore,  8 transitions to larger mean values going from  min = 20 ℎ −1 Mpc to  min = 5 ℎ −1 Mpc.The fact that these two effects are not observed in the Nseries mocks could be partially attributed to noise in the CMASS data vector, residual observational systematic effects that the mocks do not capture, or misunderstandings in the galaxy-halo connection modelling, which is mostly constrained by small scales.One of the most common systematics that can affect small-scale clustering is fiber assignment, which artificially lowers the clustering amplitude on scales smaller than the fiber collision angular scale.The Nseries mocks are already infused with fiber collisions that should closely match the BOSS data, so any important effect coming from this should also be imprinted in the multipoles measured from the mocks.
Based on the fact that the constraints from Nseries are robust against variations in  min , and that the CMASS constraints with different scale cuts are consistent to within 1, we adopt  min = 1 ℎ −1 Mpc as the baseline for our analysis.

Systematic error budget
Based on the tests described in the previous section, we use the Nseries mocks to determine the contribution to the systematic error budget coming from our emulator.This was trained on periodic boxes at a fixed redshift, but we are applying it to fit a survey that has a non-uniform footprint and radial selection, along with fiber collision effects that artificially decrease the clustering on small scales.Thus, there is the possibility of errors in addition to the emulator error term determined from the test sample when constructing the emulator.
To look for such error terms, we calculate the offset between the expected value of each cosmological parameter and the mean of the marginalized posteriors from the fit to the mean of 84 mocks, using the covariance of a 120(ℎ −3 Gpc 3 ) volume (solid contours from the right-hand side panel of Fig. 8).For the combination of DSC and the galaxy 2PCF, the offsets and their associated 2 uncertainties are: We do not find statistically significant offsets from the expected values, with the shifts for all parameters appearing well within their 2 limits.However, the level of precision to which we can carry out this test is set by the number of mocks that are available, which limits the total effective volume used in the test.We take the 2 limit of the marginalized distributions of best-fit values (left-hand side Nseries, V = 1.4 (h 3 Gpc 3 ) Nseries, V = 120 (h 3 Gpc 3 ) Figure 8. Marginalized constraints on  cdm ,  8 , and   from fits to the Nseries mock catalogues, which were calibrated onto the clustering of the CMASS NGC sample, matching its number density and geometry.The fits were performed using the baseline configuration of our analysis, consisting on the combination of the density-split and the galaxy correlation functions at scales 1.0 ℎ −1 Mpc <  < 150 ℎ −1 Mpc.Left: the dots represent the best-fit values of individual fits to 84 realizations of the Nseries mocks.The contours show the 68 and 85 per cent confidence regions of the distribution of individual fits.Right: fits to data vectors that are averaged over the 84 Nseries mocks, using a covariance rescaled to match the volume of the CMASS sample used in our main analysis (dashed-grey) or the total volume of the Nseries suite (solid-black).
of Fig. 8), divided by √ 84, as a conservative (maximum) estimate of the systematic error for each parameter, which are then added in quadrature to the associated statistical errors.We estimate this from the distribution of individual fits divided by √ 84 rather than from the fit to the mean, because the 2 limits of the latter are dominated by the emulator error on the multipoles, which is already included in the covariance when we calculate the likelihood [Eq.( 15)].In this way, we avoid double counting the emulator uncertainty, only including the limiting precision that comes from the finite number of mocks used for the test.With this setup, the systematic error added to each parameter for the combination of DSC + galaxy 2PCF combination is: These systematic errors are added in quadrature to the statistical error budget.The systematic errors are negligible compared to the statistical errors on the parameter constraints from CMASS, so the values reported in Table 3 are unaffected, up to the significant figures that are shown.

Robustness against pipeline settings
Here we explore the robustness of our clustering analysis against different choices of settings in the inference pipeline.The tests on this section are performed on the real data, and no longer using the Nseries mocks.As a reminder, our baseline configuration consists of: -Data vector: Concatenation of the monopole and quadrupole moments of the DSC auto-correlation and functions, using four quintiles (Q 0 , Q 1 , Q 3 , and Q 4 ) and the galaxy 2PCF.
-Error budget: The covariance used in the likelihood calculation includes contributions from sample variance of the data vector, and model uncertainty associated with the emulator training.
Figure 11 shows the constraints that result from varying various aspects of these settings.Starting from top to bottom, we try a uniform prior for   , finding that it shifts the means of the marginalized posteriors by less than 0.2, but degrading the constrainig power on  cdm by 15 per cent.  itself is basically unconstrained with our analysis under a uniform prior, which is the main motivation for adopting a BBN prior, following other clustering studies in the literature (e.g., Alam et al. 2017;Ivanov et al. 2020;Philcox & Ivanov 2022).
Using only the monopole of the correlation functions weakens the precision of the constraints on  cdm and  8 by 18 and 13 per cent, respectively, resulting in a 18 per cent degradation of the precision on   8 .
Using only the most underdense and overdense quintiles (Q 0 and Q 4 , respectively) degrades the precision on  cdm by 9 per cent.The precision on  8 increases by 15 per cent and its mean is shifted to slightly lower values.The increase in precision for  8 might seem counter intuitive, but can be explained by the fact that individual quintiles can sometimes lead to marginalized posteriors with slightly different mean values, which is exemplified by the rows where we fit Q 0 and Q 4 separately.This can produce wider contours when all quintiles are fitted simultaneously.Overall, the fact that the constraining power on most parameters is not severely degraded when we fit only the extreme quintiles (which effectively discards half of the DSC data set) agrees with the picture that most of the information from DSC comes from the very under-dense and over-dense regions, as suggested by Paillas et al. (2021).
DSC by itself predicts  8 = 0.768 +0.037 −0.043 , which is slightly lower than the predicted value from the 2PCF-only fit.However, both the DSC-only and 2PCF-only fits are consistent within one standard deviation with the baseline fit.
We observe that the precision on  cdm degrades by a factor of 2.4 when we allow  eff to vary, which is due to the strong correlation between these parameters.However, the precision for the other parameters remains relatively stable.If we also allow d  /d ln  and  0 to vary, it does not significantly affect the constraining power on the base ΛCDM parameters.
We find an interesting shift of the estimated  8 to lower values when we fix the assembly bias or velocity bias parameters during the fit.As shown in Appendix A, our best-fit model constrains the environment-based assembly bias parameters to be negative, meaning that galaxies preferentially form in haloes in denser environments.We observe a mild correlation between  8 and ,  vel, and  sat that is likely to explain the shifts in  8 towards lower values when removing assembly or velocity bias parameters.
Finally, we run a fit neglecting the contribution of the model uncertainty to the error budget, which can be achieved by removing C emu and C sim from the covariance used to calculate the likelihood.This has a drastic impact on the derived constraints, resulting in a precision that can be more than two times better than the baseline analysis.Although for this particular test we observe that the mean values of the recovered parameters do not shift significantly with respect to the baseline fit, we have explicitly verified that removing the emulator error results in cosmological constraints that can be biased at a more than 3 level when fitting the mean of the Nseries mocks down to 1 ℎ −1 Mpc.

DISCUSSION AND CONCLUSIONS
We have presented a clustering analysis of the DR12 BOSS CMASS galaxy sample at 0.45 <  < 0.6, using simulation-based models of the galaxy two-point correlation function (2PCF) and densitysplit clustering (DSC).Our theory framework, which is presented in detail in our companion paper (Cuesta-Lazaro et al. 2023), is based on emulators trained on high-fidelity mock galaxy catalogues, which forward model the cosmological dependence of the full shape of the galaxy 2PCF and DSC multipoles, including redshift-space and Alcock-Paczynski distortions.
It should be noted that due to the limitations of the simulation data available for training the emulator, our model fits impose a fixed prior on the acoustic scale  * measured from the CMB (Planck Collaboration et al. 2020).This is, however, a very precise and modelindependent measurement, so this prior does not significantly restrict our conclusions about the models analysed.We have validated our theory model against the Nseries mock galaxy catalogues, which were calibrated onto the clustering and selection properties of the CMASS galaxy sample, finding that we can recover unbiased cosmological constraints even using a volume of 120(ℎ −3 Gpc 3 ), which is 84 times larger than the volume examined in our data analysis.
For our base ΛCDM analysis, we find that the galaxy 2PCF constrains  cdm ,  8 , and   with a precision of 2.8, 6.1, and 3.2 per cent, respectively, using a scale range 1 ℎ −1 Mpc <  < 150 ℎ −1 Mpc.Adding the DSC multipoles using the same scale range tightens the constraining power to a precision of 1.8, 4.3 and 1.8 per cent, respectively, obtaining  cdm = 0.1201 ± 0.0022,  8 = 0.792 ± 0.034, and   = 0.970 ± 0.018.This is a factor of 1.6, 1.4, and 1.8 of improvement in precision with respect to the 2PCF-only constraints, respectively.
Combining the galaxy 2PCF and DSC multipoles, we derive   8 = 0.462 ± 0.020 at  ≈ 0.525, which is a 4.3 per cent constraint.In comparison, the main BOSS clustering analysis presented in Alam et al. (2017) derived a 8.3 per cent constraint on   8 using both Galatic caps for their 0.4 <  < 0.6 redshift bin.We obtain 1.9 times better precision using only the Northern Galactic cap and a narrower redshift bin.This improvement mainly comes from the inclusion of higher-order clustering information that is captured by DSC, and the addition of non-linear scales in the fitting.Our   8 constraint is largely consistent with Planck 2018 base-ΛCDM predictions (Planck Collaboration et al. 2020), and also agrees well with other clustering studies in the literature that use the same galaxy sample (Yuan et al. 2022b;Yu et al. 2023).Our base-ΛCDM cosmological constraints are summarized in Table 3.
We have also performed fits with single-parameter extensions to base-ΛCDM, where we vary the running of the spectral index of the primordial power spectrum (d  /d ln ), the density of massless relic neutrinos ( eff ), and the dark energy equation of state parameter  0 .Overall, we do not find compelling evidence for deviations from the base ΛCDM model, obtaining  eff = 3.02 +0.24 −0.27 , d  /d ln  = −0.005± 0.013, and  0 = −0.956+0.046 −0.041 .We have used an extended halo occupation distribution (HOD) framework to model the LRG galaxy-halo connection, using halo catalogues from the AbacusSummit suite of simulations.Our HOD constraints are largely consistent with the findings of Yuan et al. (2022a) and Yuan et al. (2022b).We constrain the minimum halo mass for hosting centrals to be log  cut = 12.65 +0.08 −0.11 , and the typical halo mass for hosting one satellite log  1 = 13.69 +0.10  −0.15 .We find signs of environment-based assembly bias, suggesting a preference for galaxies to form in haloes around denser environments.
Overall, we find that our model is able to fit the Nseries mock galaxy catalogues much more accurately than the CMASS data itself.We deem that this is not a problem specific to our emulation framework, since the mean clustering signal of the Nseries mocks, which agrees with the clustering of the AbacusSummit simulations for similar cosmologies, shows a similar offset with respect to the CMASS data.This suggests the possibility that there might be residual systematic effects in the data that are currently not captured by the mocks and are not taken into account by the weighting procedure that we adopt when calculating the clustering.
It is also worth noting that our emulator has been trained on HOD catalogues at a fixed redshift,  = 0.5, and we do not include any possible effects of redshift dependence on the modelling.Furthermore, the effective redshift of our CMASS sample,  eff = 0.525, differs slightly from the redshift of the HOD catalogues.Although the relatively narrow redshift range we are imposing in the CMASS catalogue (0.45 <  < 0.6) could alleviate some concerns regarding the redshift dependence of the theoretical model, a careful study of the impact of this simplification on the parameter constraints is going to become even more relevant when extending our framework to other datasets such as eBOSS (Dawson et al. 2016).We plan to study this in future work using the AbacusSummit lightcone simulations that have recently become available (Hadzhiyska et al. 2022).
On-going and upcoming galaxy redshift surveys, such as DESI (DESI Collaboration et al. 2016), Euclid (Laureĳs et al. 2011) and the Nancy Grace Roman Space Telescope (Green et al. 2012) will open exciting avenues to use beyond-two-point statistics for cosmology, not only in terms of improved constraints on cosmological parameters, but also in regards to our understanding of the galaxyhalo connection, the treatment of observational systematics, and the possibility of finding surprises in the data that can challenge our preconceptions about the Cosmos.
Table A1.HOD parameter constraints for the base-ΛCDM fits on the CMASS galaxy catalogue, using the combination of density-split clustering and the galaxy 2PCF, or using only the latter.In those cases where there is no constraining power other than that from the prior, we use an em dash.terms of velocity biases, we find a lower central and satellite velocity bias compared to the two previous studies.
In terms of galaxy assembly bias, we find negative  cen and  sat , which is qualitatively consistent with Yuan et al. (2022a).However, the amplitude of the inferred central assembly bias is significantly larger than that of Yuan et al. (2022a) and Yuan et al. (2022b).Physically, negative central assembly bias means that central galaxies prefer denser environments, consistent with our intuition of large red galaxies preferentially occupying highly biased cluster environments.Yuan et al. (2022b) also found a mild degeneracy between environment-based assembly bias and velocity bias.By allowing galaxies to preferentially occupy halos in denser environments, these galaxies would occupy deeper potential wells and thus have higher peculiar velocities, reducing the need to invoke additional velocity bias.This could potentially explain the lower inferred velocity bias we find.
We also notice that even though our priors on HOD parameters are fairly broad, the marginalized posteriors for some parameters hit the lower bound of the prior limits.This is the case for the velocity bias of satellites  vel,s , and the velocity bias of centrals  cen .Therefore, we advise the reader to interpret our results for these parameters as trends rather than reliable central values, error bars, or upper limits.However, our cosmology results are likely not significantly impacted, as there is weak to no correlation between these parameters and cosmological parameters.We plan to further explore different choices of priors on the HOD parameters in our upcoming work.

APPENDIX B: DEPENDENCE ON THE MODEL PARAMETERS
To explore the source of the cosmological constraining power of DSC, Fig. B1 illustrates how the model vector responds to changes in different cosmological parameters.For brevity, we only show the crosscorrelation function between the first quintile, Q 0 and the galaxy field, but we have verified that the trends for other quintiles are similar in nature.
Increasing  cdm results in a decrease in the absolute value of the amplitude of the monopole and quadrupole.There are several effects that come into play here.First, changing  cdm has an effect on the matter clustering, which affects the shape of the galaxy density PDF, and consequently the identification of density quintiles. cdm is also tightly correlated with Ω m , which determines   8 , which in turn sets the amplitude of velocity fluctuations and affects RSD.Moreover, the changes in  cdm produce AP distortions, which affect the quintilegalaxy pair separations and become more severe the more it deviates from the fiducial cosmology used to convert redshift to distances.AP distortions not only change the amplitude of the multipoles, but also produce horizontal shifts in the monopole, changing the scale of the acoustic feature.
The amplitude of the quadrupole correlates positively with  8 , which is in agreement with our expectations of linear theory regarding how the amplitude of velocity fluctuations scales with this parameter.Interestingly, the monopole amplitude remains fixed regardless of the value of  8 .This is a consequence of the way in which we populate dark matter halos with galaxies: for those HOD parameter combinations where the resulting galaxy number density is greater than the CMASS number density, ≈ 3.5×10 −4 ( ℎMpc −1 ) 3 , the catalogues are downsampled to match the target.This effectively fixes  8 in our training sample, where  is the linear galaxy bias.As we vary  8 , the galaxy bias also changes, keeping the amplitude of the monopole fixed.
Finally, we see that the spectral index of the primordial power spectrum,   , is anti-correlated with the amplitude of the monopole and quadrupole, but does not produce horizontal shifts in the profiles as in AP distortions.and the galaxy field.Each subpanel shows the monopole and quadrupole moments, rescaled by  2 to highlight the model features at large scales.The ranges that are plotted for each parameter correspond to the prior ranges used in the likelihood analysis.

Figure 1 .
Figure 1.Probability distribution function (PDF) of the galaxy overdensity measured around random query points.The overdensity field has been smoothed with a Gaussian filter of width   = 10 ℎ −1 Mpc.The colours represent the split of the PDF into quintiles.

Figure 2 .
Figure 2. Comoving number density as a function of redshift for the densitysplit quintiles and galaxies from the BOSS DR12 CMASS sample.

Figure 3 .
Figure 3. DSC and galaxy 2PCF multipoles measured from the BOSS DR12 CMASS catalogue (circles with error bars), along with the best-fit model from our emulator (solid lines with shaded bands).The columns show the multipoles of the quintile-galaxy cross-correlation function (left), the quintile autocorrelation function (middle), and the galaxy two-point correlation function (right).The lower sub-panels display the difference between the model and the data, in units of the standard deviation of the total error budget.Each colour corresponds to a different density quintile, as illustrated in Fig. 1.The 68 per cent errors of the data are estimated from 2048 realizations of the MD-Pathchy mocks and represent the expected level of sample variance for the CMASS NGC volume.The emulator uncertainty is estimated by validating the predictions against a test set of simulations with a known cosmology.
assembly bias for centrals sat U [ −1.0, 1.0]-Environment-based assembly bias for satellites following subset of simulations, which were all run with the same initial conditions: c000 PL18 ΛCDM base cosmology, matching the mean of the base_plikHM_TTTEEE_lowl_lowE_lensing likelihood.

Figure 4 .
Figure 4. Posterior probability distributions of the base-ΛCDM parameters from the fits to the BOSS DR12 CMASS clustering data.The pink contours show results from the combination of density-split clustering and the galaxy 2PCF, while the aquamarine contours show results using only the latter.We also overplot constraints taken from the Planck_TTTEEE_lowl_lowE likelihood (Planck Collaboration et al. 2020) in blue.All 2D contours show 68 and 95 per cent confidence intervals.
Figure7.Logarithm of the Bayesian evidence obtained from individual fits to the Nseries mocks (histograms) and to the BOSS DR12 CMASS data (dashed lines).We show results for two different minimum scale cuts:  min = 1 ℎ −1 Mpc (green) and  min = 50 ℎ −1 Mpc (pink).The evidence of the CMASS fit is significantly lower than any of the Nseries mocks in both cases, but this discrepancy is reduced when excluding small scales.

Figure 9 .
Figure 9. Impact of the minimum scale cut on the constraints on base-ΛCDM parameters and   8 when fitting the mean of the Nseries mocks (red circles), or the BOSS DR12 CMASS data blue squares).The horizontal dashed line shows the true cosmology of the Nseries simulations.Note that the horizontal axes do not use a linear scale.

Figure 10 .
Figure 10.Monopole of the galaxy two-point correlation function, averaged over 84 realizations of the Nseries mocks (violet circles with error bars), along with the best-fit model and its associated uncertainty (violet solid line and bands).Also shown is the monopole from the BOSS CMASS galaxy sample (grey squares with error bars).The lower sub-panel shows the difference between the Nseries mocks and the best-fit model, in units of the error bars.The dark-grey and grey shaded regions demarcate 1 and 2 offsets, respectively.

Figure 11 .
Figure 11.A comparison of the constraints on the base-ΛCDM parameters and the derived   8 value from fits run with different analysis settings.Our baseline analysis, shown at the top, was run by simultaneously fitting the galaxy 2PCF and density-split multipoles with a base-ΛCDM + extended-HOD model, using scales 1 ℎ −1 Mpc <  < 150 ℎ −1 Mpc.The other points represent variations to that baseline configuration, either by changing the data vector or the model prescription.

0 BsatFigure A1 .Figure B1 .
Figure A1.Similar to Fig. 4, but the full posterior distribution of cosmological and HOD parameters from the base-ΛCDM fits with our baseline configuration.

Table 2 .
The  2 , degrees of freedom, log-likelihood, and log-evidence for the best fits to the BOSS DR12 multipoles from the galaxy two-point correlation function (2PCF), density-split clustering (DSC) and the combination of the two.We simultaneously fit the monopole and quadrupole, using scales 1 ℎ −1 Mpc <  < 150 ℎ −1 Mpc.

Table 3 .
Parameter constraints for the base ΛCDM fits on the BOSS DR12 CMASS galaxy catalogue, using the combination of density-split clustering and the galaxy 2PCF, or using only the latter.The simple geometrical interpretation of  * means that it is one of the best-measured quantities in all cosmology (the PL18 measurement corresponding to a 0.03 per cent precision level), and this measurement is also extremely robust to changes in the cosmological model(Planck Collaboration et al. 2020).We therefore consider this to be a very well-justified prior.

Table 4 .
Constraints on single-parameter extensions to the base-ΛCDM model from the BOSS DR12 CMASS data, using the combination of densitysplit clustering and the galaxy 2PCF.
(Lange et al. 2022)2ture derived from our analysis oBautista et al. 2021;Chapman et al. 2022density-split clustering and galaxy 2PCF measurements.The base-ΛCDM prediction, based on the best-fit Planck 2018 cosmology, is shown in blue, where the darker and lighter shades represent 68 and 95 per cent confidence intervals.We also compare with other measurements from the literature at different redshifts.The data include 6dFGS(Beutler et al. 2012), eBOSS (de Mattia et al. 2021;Bautista et al. 2021;Chapman et al. 2022), as well as other studies performed on BOSS: the consensus DR12 analysis(Alam et al. 2017), a LOWZ simulation-based analysis(Lange et al. 2022), CMASS simulation-based analyses Constraints on single-parameter extensions to the base-ΛCDM model from the galaxy 2PCF + density-split clustering fits on BOSS DR12 CMASS (pink), Planck TT,TE,EE+lowl+lowE (blue), and Planck TT,TE,EE+lowl+lowE+BAO+SNe (grey).Extensions include variations in the running of the spectral index of the primordial power spectrum, d  /d ln , the effective number of relativistic species,  eff , and the dark energy equation of state parameter,  0 .
+0.046 −0.041 , which is consistent with a fiducial cosmological