# Full-Shape analysis of the power spectrum and bispectrum of DESI DR1 LRG and QSO samples

S Novell Masot, H Gil-Marín, L Verde, J Aguilar, S Ahlen, S Bailey, S Benzvi, D Bianchi, D Brooks, E Buckley-Geer, et al.

► To cite this version:

S Novell Masot, H Gil-Marín, L Verde, J Aguilar, S Ahlen, et al. Full-Shape analysis of the power spectrum and bispectrum of DESI DR1 LRG and QSO samples. Journal of Cosmology and Astroparticle Physics, 2025, 06, pp.005. 10.1088/1475-7516/2025/06/005. hal-05005592

HAL Id: hal-05005592

https://hal.science/hal-05005592v1

Submitted on 2 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Full-Shape analysis of the power spectrum and bispectrum of DESI DR1 LRG and QSO samples

To cite this article: S. Novell-Masot et al JCAP06(2025)005

View the [article online](#) for updates and enhancements.

You may also like

- Performance of the CMS high-level trigger during LHC Run 2

A. Hayrapetyan, A. Tumasyan, W. Adam et al.

- Identification of low-momentum muons in the CMS detector using multivariate techniques in proton-proton collisions at sqrt(s) = 13.6 TeV

V. Chekhovsky, A. Hayrapetyan, V. Makarenko et al.

- Muon identification using multivariate techniques in the CMS experiment in proton-proton collisions at sqrt(s) = 13 TeV

A. Hayrapetyan, A. Tumasyan, W. Adam et al.
RECEIVED: March 14, 2025

ACCEPTED: May 1, 2025

PUBLISHED: June 6, 2025

# Full-Shape analysis of the power spectrum and bispectrum of DESI DR1 LRG and QSO samples

**S. Novell-Masot et al.**

Full author list at the end of the paper

E-mail: sergi.novell@icc.ub.edu, hectorgil@icc.ub.edu, liciaverde@icc.ub.edu

**ABSTRACT:** We present the first joint analysis of the power spectrum and bispectrum using the Data Release 1 (DR1) of the Dark Energy Spectroscopic Instrument (DESI), focusing on Luminous Red Galaxies (LRGs) and quasars (QSOs) across a redshift range of $0.4 \le z \le 2.1$. By combining the two- and three-point statistics, we are able to partially break the degeneracy between the logarithmic growth rate, $f(z)$, and the amplitude of dark matter fluctuations, $\sigma_{s8}(z)$, which cannot be measured separately in analyses that only involve the power spectrum. In comparison with the (fiducial) Planck ΛCDM cosmology we obtain $f/f^{\text{fid}} = \{0.888_{-0.089}^{+0.186}, 0.977_{-0.220}^{+0.182}, 1.030_{-0.085}^{+0.368}\}$, $\sigma_{s8}/\sigma_{s8}^{\text{fid}} = \{1.224_{-0.133}^{+0.091}, 1.071_{-0.163}^{+0.278}, 1.000_{-0.223}^{+0.088}\}$ respectively for the three LRG redshift bins, corresponding to a cumulative 10.1% constraint on $f$, and of 8.4% on $\sigma_{s8}$, including the systematic error budget. Additionally, we obtain constraints for the ShapeFit compressed parameters describing the isotropic scaling parameter, $\alpha_{\text{iso}}(z)$, the Alcock-Paczyński parameter, $\alpha_{\text{AP}}(z)$, the combined growth of structure parameter $f\sigma_{s8}(z)$, and the combined shape parameter $m(z) + n(z)$. Their cumulative constraints from our joint power spectrum-bispectrum analysis are respectively $\sigma_{\alpha_{\text{iso}}} = 0.9\%$ (9% improvement with respect to our power spectrum-only analysis); $\sigma_{\alpha_{\text{AP}}} = 2.3\%$ (no improvement with respect to power spectrum-only analysis, which is expected given that the bispectrum monopole has no significant anisotropic signal); $\sigma_{f\sigma_{s8}} = 5.1\%$ (9% improvement); $\sigma_{m+n} = 2.3\%$ (11% improvement). These results are fully consistent with the main DESI power spectrum analysis, demonstrating the robustness of the DESI cosmological constraints, and compatible with Planck ΛCDM cosmology.

**KEYWORDS:** cosmological parameters from LSS, power spectrum, redshift surveys

ARXIV ePRINT: 2503.09714

<figure><img src="image_2.png" /></figure>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Overview of the DESI DR1 catalogues and mocks</td><td>3</td></tr><tr><td>2.1</td><td>Mocks</td><td>6</td></tr><tr><td>3</td><td>Baseline analysis</td><td>7</td></tr><tr><td>3.1</td><td>Standard compression and ShapeFit</td><td>9</td></tr><tr><td>3.2</td><td>Power spectrum modelling</td><td>10</td></tr><tr><td>3.3</td><td>Bispectrum modelling</td><td>12</td></tr><tr><td>3.4</td><td>Likelihood and covariance estimation</td><td>15</td></tr><tr><td>4</td><td>Systematic error budget</td><td>15</td></tr><tr><td>4.1</td><td>Theoretical systematics</td><td>16</td></tr><tr><td>4.2</td><td>Fibre assignment and imaging systematics</td><td>18</td></tr><tr><td>4.3</td><td>HOD systematics</td><td>21</td></tr><tr><td>4.4</td><td>Systematics due to fiducial cosmology assumption</td><td>24</td></tr><tr><td>5</td><td>Results</td><td>26</td></tr><tr><td>5.1</td><td>Unblinded data</td><td>28</td></tr><tr><td>5.2</td><td>Statistical error validation</td><td>31</td></tr><tr><td>6</td><td>Conclusions</td><td>34</td></tr><tr><td>A</td><td>Perturbation theory modelling of the power spectrum and bispectrum</td><td>36</td></tr><tr><td>B</td><td>Blinded data</td><td>38</td></tr><tr><td>C</td><td>Nuisance parameters</td><td>40</td></tr><tr><td>D</td><td>Covariances and additional tables</td><td>41</td></tr><tr><td></td><td>Author List</td><td>52</td></tr></table>

## 1 Introduction

Cosmology is in its precision era, with increasingly larger datasets providing us with huge statistical power to test our cosmological models of the Universe. The Λ-Cold Dark Matter (ΛCDM) model has proven remarkably successful in explaining observations across a broad range of epochs. However, fundamental questions persist about the nature of the dark component of the universe (dark matter and dark energy), which constitutes the majority of the Universe's energy density. In addition, several tensions, such as the Hubble tension [1], persist between early and late Universe observations, thus suggesting that further investigation of the possible inconsistencies of ΛCDM is necessary.
Galaxy redshift surveys have been instrumental in probing the large-scale structure (LSS) of the Universe. Surveys such as the Sloan Digital Sky Survey (SDSS) [2], the 2-degree Field Galaxy Redshift Survey (2dFGRS) [3], the Baryon Oscillation Spectroscopic Survey (BOSS) [4] and extended BOSS (eBOSS) [5] have allowed the cosmology community to, sequentially, lay the foundations, develop the methodology, and obtain precise constraints on the baryon acoustic oscillation (BAO) feature and the growth of structure.

The Dark Energy Spectroscopic Instrument (DESI) [6–21] builds on this legacy, targeting over 40 million objects, including galaxies, quasars, and Lyα absorbers, thus mapping the cosmic web across a wide range of redshifts ($z \sim 0$ to $z \sim 4$). This dataset spans a volume several times larger than previous surveys with enough signal to go beyond two-point statistics, the primary workhorse of cosmology analyses to date. During the course of five years of observations, DESI will cover an effective volume of $V_{\text{eff}} \sim 50–60 \, \text{Gpc}$.

The bispectrum, which is the natural next-order statistic, quantifies the correlations among triplets of points (in an analogous way as the power spectrum is obtained by correlations of pairs). Given that it is a higher-order statistic, the bispectrum encloses information about the non-Gaussian features of the matter distribution, probing non-linear structure formation, non-linear galaxy bias, and also potential signatures of primordial non-Gaussianities [22–24]. Hence the bispectrum complements the power spectrum by helping to break degeneracies among cosmological parameters [25–42], as will be discussed in what follows.

Application to real data is however challenging and computationally intensive: the bispectrum data-vector is large, with correlations between elements, its signal-to-noise is low, the covariance matrix is non-diagonal. Additionally, the modelling of the signal and the potential systematics is highly complex. This is why, compared to the extensive theoretical efforts, application of the bispectrum to real data to constrain cosmology has been somewhat limited to date, with few attempts by very few research groups [41, 43–54].

In this work, we perform the first joint power spectrum-bispectrum analysis using the DESI Data Release 1 (DR1) [55], for both the Luminous Red Galaxy (LRG) and the quasar (QSO) samples. These samples span across the redshift range $0.4 \le z \le 2.1$, thus probing both early and intermediate epochs of structure formation. We model the bispectrum signal employing the GEO-FPT model, which has been calibrated and validated using a broad range of simulations [56], while our power spectrum analysis is based on the renormalized perturbation theory model (RPT, [57]) and resembles the one employed in the BOSS and eBOSS analyses [58, 59]. To accurately quantify both our statistical and systematic errors, we employ sets of realistic mock catalogues that allow us to estimate the covariance matrix and quantify observational and theoretical systematics.

We derive constraints on cosmological parameters within the template-based *ShapeFit* approach of [60]. This approach is referred to as a type of Full-Shape analyses in the DESI papers [61], as it exploits information on the shape of the power spectrum (and bispectrum). Other Full-Shape approaches, such as the Full Modelling (also known as direct fit) employed in other DESI papers are not considered in this work. Within ShapeFit, the *shape* parameters *m* and *n*, enclosing information about the shape of the matter power spectrum, are added to the standard parameters of interest in template-based methods: the dilation parameters along and across the line-of-sight, $\alpha_\|$ and $\alpha_\perp$, the logarithmic growth rate of structure, *f*, and the amplitude of matter fluctuations, $\sigma_{s8}$.
As detailed in [60], we adopt $\sigma_{s8}$ (rather than the more common parameter $\sigma_8$) for parametrising the smoothed amplitude of matter fluctuations. The essential difference is that $\sigma_8$ smooths the fluctuations on the fixed scale of $8 h^{-1}$Mpc (thus $\sigma_8$ varies with the change of scales via $\alpha_\parallel, \alpha_\perp$) while $\sigma_{s8}$ defines the smoothing scale as to be independent of changes in scale because it takes a fixed smoothing scale at the template cosmology, $s \equiv r_d/r_d^{\text{fid}}$:

$$ \sigma_{s8} \equiv \sigma(R = s \cdot 8 h^{-1} \text{Mpc}), \qquad (1.1) $$

with $s$ depending on the cosmology at each step of the fitting process. This change in definition does not alter the fitting part of the analysis, only the interpretation of results. In the case where the best-fit dilation parameters $\alpha_\parallel, \alpha_\perp$ are equal to 1 and the true sound horizon scale is equal to the fiducial sound horizon scale, $\sigma_{s8}$ coincides with $\sigma_8$.

As already noted in [60], the shape parameters $m$ and $n$ are very degenerate and we are unable to constrain them independently, thus we only constrain their combination $m+n$. Analogously, in analyses only involving two-point statistics, the parameters $f$ and $\sigma_{s8}$ have a very strong degeneracy, which is why it is common to only consider their product, $f\sigma_{s8}$. However, the inclusion of the bispectrum breaks the degeneracy<fnref n="1" /> and $f$ and $\sigma_{s8}$ can be inferred separately [52, 56].

The structure of this paper is as follows. In section 2 we review the main aspects of the DESI DR1 data and mock catalogues. In section 3 we describe our adopted modelling and baseline choice for the power spectrum and bispectrum analysis and compare it with the standard approach from the DESI collaboration. In section 4 we quantify the impact of the main sources of systematic errors. In section 5 we present our results, which we also compare with the official DESI full-shape power spectrum analysis of [61]. We conclude in section 6.

## 2 Overview of the DESI DR1 catalogues and mocks

The DESI DR1 galaxy data [55] encompasses spectroscopic observations made by the DESI telescope in the time period between May 14, 2021 and June 14, 2022. One of the main improvements of the DESI telescope with respect to previous experiments is its structure into ten ‘petals’, that in total comprise 5000 fibres, each guided by a robotic positioner. According to observing conditions, each fibre gets assigned a target, whose light is redirected to one of the ten spectrographs. This improvement has resulted in that the first year of observations (of total effective volume 19.5 Gpc³) has a statistical power which surpasses that of two decades of the Sloan telescope observations at Apache Point [62].

There are four main observed objects: Bright Galaxy Survey (BGS) [63], Luminous Red Galaxies (LRG, divided into three redshift bins LRG1, LRG2, LRG3) [14], Emission Line Galaxies (ELG) [64] and quasars (QSO) [65]. In all cases, the observations are divided into the North Galactic Cap (NGC) and the South Galactic Cap (SGC). Since both regions have been tested to be compatible, as expected, we consider the joint NGC+SGC data in all cases in this paper, as is also done in the DESI main analyses [17, 20, 61, 66]. We use the three LRG redshift bins, which have the most signal-to-noise [14], and the QSOs.<fnref n="2" /> The ELG sample

<footnote n="1">¹This is because $f$ and the bias parameters $b_1, b_2$, all degenerate with $\sigma_{s8}$ in the power spectrum, are related differently in the power spectrum than in the bispectrum, as is seen in appendix A.</footnote>

<footnote n="2">²We do not include the BGS sample due to its lower signal.</footnote>
<table><thead><tr><th>Tracer</th><th>Redshift range</th><th>N<sub>tracer</sub></th><th>z<sub>eff</sub></th><th>P<sub>0</sub> [(Mpc h<sup>-1</sup>)<sup>3</sup>]</th><th>V<sub>eff</sub> [Gpc<sup>3</sup>]</th><th>Used in this work</th></tr></thead><tbody><tr><td>BGS</td><td>0.1 – 0.4</td><td>300,017</td><td>0.295</td><td>~ 9.2 × 10<sup>3</sup></td><td>1.7</td><td>No</td></tr><tr><td>LRG1</td><td>0.4 – 0.6</td><td>506,905</td><td>0.510</td><td>~ 8.9 × 10<sup>3</sup></td><td>2.6</td><td>Yes</td></tr><tr><td>LRG2</td><td>0.6 – 0.8</td><td>771,875</td><td>0.706</td><td>~ 8.9 × 10<sup>3</sup></td><td>4.0</td><td>Yes</td></tr><tr><td>LRG3</td><td>0.8 – 1.1</td><td>859,824</td><td>0.920</td><td>~ 8.4 × 10<sup>3</sup></td><td>5.0</td><td>Yes</td></tr><tr><td>ELG1</td><td>0.8 – 1.1</td><td>1,016,340</td><td>0.955</td><td>~ 2.6 × 10<sup>3</sup></td><td>2.0</td><td>No</td></tr><tr><td>ELG2</td><td>1.1 – 1.6</td><td>1,415,687</td><td>1.317</td><td>~ 2.9 × 10<sup>3</sup></td><td>2.7</td><td>No</td></tr><tr><td>QSO</td><td>0.8 – 2.1</td><td>856,652</td><td>1.491</td><td>~ 5.0 × 10<sup>3</sup></td><td>1.5</td><td>Yes</td></tr></tbody></table>

**Table 1.** Summary of the DESI DR1 tracers main properties, including redshift ranges, tracer counts, effective redshift $z_{\text{eff}}$, power spectrum monopole amplitude at the reference scale of $k = 0.14 \text{ hMpc}^{-1}$, $P_0$, and effective volume $V_{\text{eff}}$. The effective redshift and volume, $z_{\text{eff}}$, $V_{\text{eff}}$ are defined in equation (2.1). In this work, we do not consider the BGS nor the ELG samples, focusing only on the three LRG redshift bins together with the QSO redshift bin.

features large imaging systematics [61, 67], which are not easily treated within our pipeline, so we do not use these samples in this work. The details of the different samples are shown in table 1. The mathematical definitions used for the quoted quantities $z_{\text{eff}}$, $V_{\text{eff}}$ in table 1 are

$$z_{\text{eff}} = \frac{\sum_{i>j} w_i w_j (z_i + z_j)/2}{\sum_{i>j} w_i w_j} \quad [59]; \quad V_{\text{eff}} = \int \left[ \frac{\bar{n}(\mathbf{r}) P(k)}{1 + \bar{n}(\mathbf{r}) P(k)} \right]^2 d^3 r \quad [68], \quad (2.1)$$

where, respectively $w_i, w_j$ are the galaxy weights, and $\bar{n}$ is the density of the sample.

In order to account for the effect of the survey geometry, a set of random catalogues are produced, consisting of a Poissonian distribution covering the same survey footprint as the data (both angular and radial), but with no structure. As it is standard practice, these random catalogues are combined with the data, following the FKP estimator approach [69]. For both data and mocks, the power spectra and bispectra measurements are obtained with the RUSTICO code<fnref n="3" /> [59], taking into account both the observational systematics weights and the so-called FKP weights [69].<fnref n="4" /> Figure 1 shows the measurements of the power spectrum and bispectrum multipoles of our baseline analysis (for details and motivation see appendix B), together with the ratio between the model best fit and the measurements.

Additionally, we adopt the blinding policy of the DESI key projects, according to which we only analyse the data once all analysis choices (e.g., modelling, range of scales, choice of tracers) are set and ‘frozen’. Until then, we only consider the mocks of the tracers of interest — presented in section 2.1 below — and the blinded data. The blinding algorithm [70] has two components, which shift the position of the BAO peak (and thus the recovered dilation parameters $\alpha_\parallel, \alpha_\perp$) and the growth rate parameter $f$ respectively. Additionally, there is a third component of the blinding in DESI, affecting the signature of the local primordial non-Gaussianity parametrised by $f_{\text{nl}}$ [71]. In the context of our present analysis, where we kept fixed the $f_{\text{nl}}$ term to 0 in our modelling, the $f_{\text{nl}}$ blinding that is part of the blinded DR1

<footnote n="3">³https://github.com/hectorgil/Rustico.</footnote>

<footnote n="4">⁴The FKP weights ensure the minimum variance of the power spectrum at a certain chosen scale of reference, $0.14 \text{ hMpc}^{-1}$, by balancing contributions between regions with different number densities.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Power spectrum (top section) and bispectrum (bottom section) measurements and comparison with the respective model best fit. Each column corresponds to a tracer. First row: In each panel, power spectrum monopole (top line), quadrupole (middle) and hexadecapole (bottom). The points and the shaded region are the measurements with their associated error, while the solid lines show the best fit of our joint $P+B$ analysis (section 5.1). The hexadecapole of the QSO is not used in this analysis as discussed in appendix B); Middle three rows: offset between the best fit of our $P+B$ analysis and the respective measurement, divided by the error bar of the measurement (shaded areas enclose the $2.5\sigma$ region). Bottom panels: bispectrum monopole, and using the same convention as above.</figcaption></figure>
catalogues causes a shift in the ShapeFit $m+n$ parameter as these two parameters are very degenerate [72]. Refs. [70, 73] validated the blinding scheme, demonstrating how the recovered BAO and redshift-space distortions (RSD) parameters are shifted coherently with the expected blinding shift. This study was also extended for the power spectrum and bispectrum RSD analysis in [74], validating this catalogue blinding scheme also for bispectrum analyses.

## 2.1 Mocks

We have two main needs for using mock catalogues: 1) to assess and quantify the importance of the various systematic errors in the analysis, and 2) to estimate the covariance matrix of the data. Since there is a trade-off between accuracy and computational cost in producing simulations, two distinct sets of mocks have been generated within the DESI collaboration to address these two separate needs.

- **AbacusSummit.** A suite of high-resolution N-body gravity simulations produced with the ABACUS N-body code [75], each containing 6912³ in a physical volume of $V_{\text{box}} = 8 \, (\text{Gpc}/h)^3$. Although there are a total of 97 cosmologies in the AbacusSummit suite, we only employ the 25 so-called 'base' realisations in all tests. The halos are identified using the COMPASSO halo finder, and they are populated with the appropriate tracers according to a halo occupation distribution (HOD) [76, 77]. These mocks, which we hereafter refer to as Abacus, represent the expected clustering of DESI tracers with high fidelity. Hence, we use them to perform tests in order to obtain the systematic error budget (see section 4) and specify the analysis set-up choices (see section 3).

- **EZmocks.** These approximate mocks, which need to be calibrated to the Abacus mocks, directly generate the density field using the Zel'dovich approximation instead of performing the full N-body evolution. After that, the density field is populated by the corresponding tracers using an effective bias model. Since they can be obtained at a much reduced computational cost, a set of 1000 such realisations of volume $V_{\text{box}} = 6^3 \, (\text{Gpc}/h)^3$ is produced for each galaxy sample, which are used to estimate covariance matrices. As explained in detail in [61] and seen in figure 2, the EZmocks do not fully reproduce the statistics of the DESI DR1 data, which, if unaccounted for, results in an underestimate of the errors. To correct for this, as prescribed by [61], we rescale the covariance matrix of each redshift sample by the factors in table 3 of [61].

In both cases, there are two main versions of the mocks, usually labelled as respectively first and second generation mocks. We use the second generation mocks, whose HOD was calibrated with the full DESI Early Data Release (EDR) sample [16].

The cosmology of both sets of mocks corresponds to the Planck 2018 best fit, which we may denote as 'base' or c000 (following the nomenclature from [75]) throughout this paper: the densities of cold dark matter and baryonic matter being respectively $\Omega_{\text{cdm}}h^2 = 0.12$ and $\Omega_b h^2 = 0.02237$; the amplitude of the dark matter fluctuations at $z=0$ and the slope of the primordial power spectrum take the values of $\sigma_{\text{s8}} = 0.81135$ and $n_s = 0.9649$, respectively; the dimensionless Hubble parameter is $h = 0.6736$, and the dark energy equation of state corresponds to a cosmological constant, so $w_0 = -1$ and $w_a = 0$. This is also the fiducial cosmology in our analysis.
Then, the periodic boxes (both Abacus and EZmocks) are transformed into samples that account for the observed survey geometry, which we refer to as cutsky samples. In this process the Cartesian coordinates are transformed first into sky coordinates (redshift, right ascension and declination) and then the tracers are filtered so as to obtain a density matching the DESI footprint and radial selection;<fnref n="5" /> finally, the process of fibre assignment of the DESI instrument is simulated by running the FIBREASSIGN pipeline. This last step reproduces the effect of missing galaxies during the DESI observing strategy, which is caused by the limitation of placing the robotic fibre positioners on two objects which are angularly very close [19].

Depending on how the fibre assignment is simulated, three types of mocks are produced, which serve various purposes:

- **Complete.** Without fibre assignment incompleteness, so these can be considered 'uncontaminated' mocks.

- **AltMTL [78]**. The fibre assignment is implemented in the mocks using the same strategy as in the observed data. Although it is the most realistic approach, the computational cost associated only allows the AltMTL pipeline to be used with the sets of 25 Abacus mocks.

- **Fast-fibreassign (FFA).** This approach approximates the fibre assignment procedure, using averaged targeting probabilities from the data. Given that this procedure is significantly faster than the AltMTL pipeline, the FFA code is used to obtain the sets of 1000 EZmocks with fibre assignment effects.

For both AltMTL and FFA mocks, the clustering and the associated randoms are run through the corresponding fibre assignment pipeline. We plot in figure 2 the comparison between the mean data-vectors of the Abacus AltMTL, Complete, and EZmocks FFA for the LRG2 bin, together with the corresponding DESI DR1 measurement.

As we will show in section 4, we employ the complete mocks to assess the error associated with the modelling, while we can quantify the systematics originated by the fibre assignment process from the difference between the parameters recovered from the complete vs AltMTL mocks.

## 3 Baseline analysis

In this work we opt to perform a full-shape analysis of the DESI DR1 power spectrum and bispectrum data using the compression technique of ShapeFit. This approach allows us to compress all the information contained in the power spectrum and bispectrum, in a given redshift bin, into a set of parameters that do not depend on the specific choice of cosmological model (ΛCDM, wCDM, kCDM, etc.). In the companion paper [79] these compressed parameters are interpreted in the light of multiple models, bypassing the lengthy

<footnote n="5">⁵This leads to the radial integral constraint effect, described and studied in [19] for two-point statistics. As seen in [39, 56], most of the added information from the bispectrum statistic comes from linear to mildly non-linear scales, so small shifts on scales of $k \sim 0.02 h \text{Mpc}^{-1}$ produce a negligible effect in our analysis. This is also seen in figure 7, where the mode correction term does not produce a significant change in the joint power spectrum-bispectrum cosmological parameters posteriors.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Power spectrum and bispectrum as measured from the LRG2 main mocks (EZmocks FFA, Abacus AltMTL and complete), together with the DESI DR1 measurement from the same sample. We show the three power spectrum multipoles in the upper panel and the bispectrum monopole in the lower panel. For the mocks, we plot the mean of each data-vector across all available realisations (25 for Abacus and 1000 for EZmocks). Note that the normalisation of each of these displayed power spectra and bispectra is an arbitrary quantity that is only required (for unbiased measurements) to be consistent between each of these datasets and their corresponding window matrix. This means that across samples each normalization is in principle arbitrary and unrelated. For an insightful comparison, in this figure we (arbitrarily) choose to display *all* power spectra and bispectra normalization relative to the one of the AltMTL mocks.</figcaption></figure>

process of refitting again the power spectrum and bispectrum. This approach is complementary to the direct fit approach where the parameters of a given model are directly fitted to the summary statistics. We refer the reader to section 4.2 of [61] (and references therein) for a detailed description of these two methodologies.

Our baseline choice for the analysis is the combination of the power spectrum multipoles and the bispectrum monopole of LRGs and QSO samples as follows. For the LRG samples we consider the first three non-zero multipoles: monopole, quadrupole and hexadecapole in combination with the bispectrum monopole, $P_{024} + B_0$; whereas for the QSO sample we only consider the power spectrum monopole and quadrupole in combination with the bispectrum monopole, $P_{02} + B_0$. This choice is guided by the rationale of minimising the statistical errors while keeping the systematics (as quantified on mocks) under control, see appendix B for more details.
## 3.1 Standard compression and ShapeFit

As described in [60, 80, 81], the ShapeFit compression approach increases the information extracted from summary statistics with respect to the ‘standard compression’ technique used previously in BOSS and eBOSS analyses [58, 59, 82–89].

The standard compression approach encapsulates the information content of the data into three main parameters: the dilation parameters along and across the line of sight, $\alpha_\parallel$ and $\alpha_\perp$, respectively; and $f\sigma_{s8}$, the product of the (smoothed) amplitude of dark matter fluctuations ($\sigma_{s8}$) and the logarithmic growth rate ($f$). The parameters $\alpha_\parallel, \alpha_\perp$ describe the distance dilation distortions caused by the use of a fixed fiducial cosmology when constructing the catalogue, and the isotropic shift of the BAO peak through the value of the sound horizon scale on the fixed-template fiducial cosmology:

$$ \alpha_\parallel(z) = \frac{D_H(z)r_d^{\text{fid}}}{D_H^{\text{fid}}(z)r_d}; \quad \alpha_\perp(z) = \frac{D_M(z)r_d^{\text{fid}}}{D_M^{\text{fid}}(z)r_d}, \qquad (3.1) $$

where $D_H$ is the Hubble distance, $D_H(z) \equiv c/H(z)$, where $H(z)$ is the Hubble expansion history, and $D_M$ is the comoving angular diameter distance. Here $r_d$ denotes the sound horizon scale at radiation drag and the superscript ‘fid’ stands for the fiducial cosmology chosen for creating the catalogue (through the redshift-to-distance transformation) and for the fixed-template cosmology. In section 5 the main results will be presented in terms of an alternative re-parametrization of these dilation parameters,

$$ \alpha_{\text{iso}}(z) \equiv [\alpha_{\parallel}(z)\alpha_{\perp}(z)^2]^{1/3} = \frac{D_V(z)r_d^{\text{fid}}}{D_V^{\text{fid}}(z)r_d}, \quad \alpha_{\text{AP}}(z) \equiv \frac{\alpha_{\parallel}(z)}{\alpha_{\perp}(z)} = \frac{D_H(z)D_M^{\text{fid}}(z)}{D_H^{\text{fid}}(z)D_M(z)}, \qquad (3.2) $$

where $D_V(z) \equiv [zD_H(z)D_M(z)^2]^{1/3}$ is sometimes referred as the isotropic BAO distance. In this form, the effect of the absolute size of the BAO peak position in the fixed-template approach, $\propto r_d/r_d^{\text{fid}}$, which only enters into the isotropic dilation scale, is separate from a purely anisotropic component, $\alpha_{\text{AP}}$, known as the Alcock-Paczyński parameter [90].<fnref n="6" />

The $f\sigma_{s8}$ parameter quantifies the growth of structure and RSD and it is the natural parameter combination constrained by the two-point statistics: there is a strong degeneracy between $f$ and $\sigma_{s8}$. Hence, when only the power spectrum multipoles are employed, we fit for $f$ with $\sigma_{s8}$ fixed at its fiducial value, $\sigma_{s8}^{\text{fid}}$, and later re-interpret the fit of $f$ as $f\sigma_{s8}$.

Furthermore, adding a three-point statistic such as the bispectrum might allow the $f\sigma_{s8}$ degeneracy to be lifted. Therefore, in this work, cosmological constraints involving the bispectrum consider $f$ and $\sigma_{s8}$ separately in all the LRG redshift bins. The QSO sample, however, has a weak signal and thus we instead use the joint parameter $f\sigma_{s8}$ for the QSO tracers.

Additionally to the set of parameters $\{\alpha_{\text{iso}}, \alpha_{\text{AP}}, f, \sigma_{s8}\}$, the ShapeFit approach used in both this work and in the full-shape DESI power spectrum analysis [60, 61], includes a further parameter governing the combined shape parameter: $m+n$. The purpose of $m+n$ is to capture information on the matter-radiation equality present in the power spectrum

<footnote n="6">⁶Note however, that even when considering the isotropic dilation parameter at different redshifts, the ratio between $\alpha_{\text{iso}}$ at different redshifts is independent of $r_d/r_d^{\text{fid}}$.</footnote>
(through $m$), but also to account for variations of the primordial spectral index (through $n$). It does so by modifying the linear power spectrum $P_L$ as,

$$P_L(k) \rightarrow P_L(k) \exp \left\{ \frac{m}{a} \tanh \left[ a \ln \left( \frac{k}{k_p} \right) \right] + n \ln \left( \frac{k}{k_p} \right) \right\}, \quad (3.3)$$

with $k_p$ and $a$ set as $k_p = \pi/r_d^{\text{fid}}; a = 0.6$, as justified in [60]. Due to the strong degeneracy between $m$ and $n$ [60, 91], we opt to define the variables $m' \equiv m|_{n=0}, n' \equiv n|_{m=0}$, which satisfies $m' \approx n' \approx m+n$, to which we refer as the combined shape parameter. For this reason we fix $m$, and only vary $n$. In the interpretation step, this $n$-at-a-fixed-$m$ parameter is equivalent to the combination $m+n$, to which we refer as the combined shape parameter. This re-interpretation of the fit can also be seen as the linear re-parametrization $\{m, n\} \rightarrow \{m+n, m-n\}$, where $m+n$ is well-constrained, whereas $m-n$ is a poorly constrained parameter that we decide to fix to a fiducial choice $m^{\text{fid}} - n^{\text{fid}} \equiv 0$.

Indeed, we have checked that performing parameter inference by fitting either $n'$ or $m'$ produces indistinguishable results.

The variation of $m+n$ in the fit modifies the $\sigma_{s8}$ parameter such that

$$\sigma_{s8} = \sigma_{s8}^{\text{fid}} A^{1/2} \exp \left( \frac{m+n}{2a} \tanh \left[ a \ln \left( \frac{r_d^{\text{fid}}}{8h^{-1}\text{Mpc}} \right) \right] \right), \quad (3.4)$$

with $\sigma_{s8}^{\text{fid}} \equiv \sigma_8^{\text{fid}}$, and $A = A_{\text{sp}}/A_{\text{sp}}^{\text{fid}}$, where

$$A_{\text{sp}} = \left( \frac{r_{\text{d}}^{\text{fid}}}{r_{\text{d}}} \right)^3 P_{\text{L, nw}} \left( k_p \frac{r_{\text{d}}^{\text{fid}}}{r_{\text{d}}} \right), \quad (3.5)$$

with $P_{\text{L, nw}}$ being the broadband (without BAO) linear dark matter power spectrum.

In addition to the cosmological parameters of interest, $\{\alpha_{\text{iso}}, \alpha_{\text{AP}}, f, \sigma_{s8}, m+n\}$, we consider the following nuisance parameters: the first and second-order bias parameters $b_1$ and $b_2$; the deviations from Poissonian shot-noise for the power spectrum and bispectrum, respectively $A_{\text{P}}, A_{\text{B}}$ [45, 59]; and the Fingers-of-God damping factors, also varied independently for the power spectrum ($\sigma_{\text{P}}$) and bispectrum ($\sigma_{\text{B}}$). The remaining non-linear bias parameters of the two-loop bias expansion are fixed according to the hypothesis of local Lagrangian bias: the tidal bias $b_{s2}$ and third-order bias $b_{3\text{nl}}$ are set as functions of $b_1$, respectively $b_{s2} = -\frac{4}{7}(b_1 - 1)$ and $b_{3\text{nl}} = 32/315(b_1 - 1)$ [80, 92, 93]. This choice of bias expansion, compatible with the official analysis of BOSS and eBOSS [58, 59], is different than the one used in the main DESI collaboration for full-shape power spectrum. The main DESI collaboration analysis [61] includes counter-terms and stochastic parameters through the effective field theory formalism (EFTofLSS), which modify the impact of $b_{3\text{nl}}$ and $b_{s2}$ on the shape of the power spectrum.

We comment on the details of our modelling and analysis choices of the power spectrum and bispectrum in the following two subsections.

## 3.2 Power spectrum modelling

We model the power spectrum, (both for matter, $P_{\text{NL}}$, and for galaxies, $P_g$) at two-loop in renormalized perturbation theory (RPT),<fnref n="7" /> following [45, 94], which closely resembles

<footnote n="7">⁷The implementation of this code can be found in https://github.com/hectorgil/PTcool.</footnote>
the approach taken by the eBOSS collaboration to model LRGs in Fourier space [59]. In appendix A we provide a brief review, referencing the original sources of this formalism. This model, based on standard perturbation theory (SPT hereafter), differs from the main modelling choices of the DESI collaboration [61, 91, 95–98], which are based on the EFTofLSS framework. The main effect of our choice is that the power spectrum can be used for scales not smaller than $k = 0.15 \, [h\,Mpc^{-1}]$ in the redshifts corresponding to the LRG tracers, while the EFT approaches are used up to $k = 0.20 \, [h\,Mpc^{-1}]$. However, it has been tested (figure 20 in appendix B of [61] for the DR1 Data; and figure 4 of [96] for the ABACUS N-body mocks) that the information gain of going beyond $k = 0.16 \, [h\,Mpc^{-1}]$ is small. As we show later in figure 3, the difference between using RPT or EFT models (FOLPS in that case) is very small even when the mock volume is of the order of 200 [$Gpc^{-1}$]$^3$. This result is in agreement with previous results such as those displayed in [80].

We consider in all cases two different sets of power spectrum multipoles: the monopole and quadrupole, $P_{02} = \{P_0, P_2\}$ (as in the official DESI analysis [61]), and the additional inclusion of the hexadecapole, $P_{024} = \{P_0, P_2, P_4\}$. In appendix B we test the effect of adding the hexadecapole, both on the constraints and on the systematic error budget. Based on those findings, we decide (before unblinding the data) to include the power spectrum hexadecapole as our baseline model in the three LRG bins,<fnref n="8" /> but not in the QSO bin (see appendix B).

The effect of the survey window function is modelled as in [59, 82, 99]. This approach, which follows that of the BOSS and eBOSS collaborations, has some differences with the DESI official pipeline: the DESI collaboration implements a cut in the contribution of galaxy pairs separated by less than 0.05 deg (in a procedure labelled as $\theta$-cut), and a rotation of the data-vector, window matrix and covariance to better handle the resulting mode coupling [61]. The $\theta$-cut method [100] mitigates fibre assignment systematic errors (which will be described in section 4.2), and is increasingly necessary when small scales are included in the analysis, as in the official DESI full-shape analysis, which has $0.02 < k \, [h\,Mpc^{-1}] < 0.20$. In our analysis we consider larger scales ($k < 0.15 \, [h\,Mpc^{-1}]$) for the power spectrum and $k < 0.12 \, [h\,Mpc^{-1}]$ for the bispectrum), so we opt to not use the $\theta$-cut correction. We quantify the residual effects of unmitigated fibre assignment systematics in section 4.2.

### 3.2.1 Comparison of the performance of power spectrum models: RPT vs EFT

For completeness, we display here a direct comparison of the recovered parameters, from the power spectrum monopole and quadrupole only, for the two-loop RPT model used in this paper and Folps, which is an implementation of one of the EFTofLSS codes used in the DESI full-shape analysis [61, 91]. The left panel of figure 3 displays the posteriors of the ShapeFit parameters obtained by performing a power spectrum monopole and quadrupole fit to the

<footnote n="8">⁸The power spectrum hexadecapole has not been included in the official DESI power spectrum analysis due to its contribution to projection effects, and to the fact that within ΛCDM the information gain provided by $P_4$ is small. In the case of our power spectrum model, we proceed as in its previous application to BOSS and eBOSS [58, 59] and include the hexadecapole for the LRG tracers, which does not worsen the projection effects nor the systematic error with respect to the analysis without $P_4$, as seen in appendix B.</footnote>
mean of the 25 Abacus LRG cubic mocks at $z = 0.8$,<fnref n="9" /> with the covariance rescaled to the total of the 25 mocks (volume of 200 [Gpc h⁻¹]³). In this comparison, we do not employ the hexadecapole signal to match the setup baseline choice of the main DESI full-shape analysis.

The figure displays two setup variations of Folps, one with the so-called ‘maximal freedom’ (empty contours in light orange), and ‘minimal freedom’ (filled contours in red). As described in [91], the maximal freedom leaves the bias parameters $b_{s2}$ and $b_{3nl}$ free, and the minimal freedom sets them to the Local Lagrangian predictions (or co-evolution relations, see eq. (A.2)). In both cases, the range of fitted scales is $0.02 < k [hMpc^{-1}] < 0.18$.<fnref n="10" /> Our RPT model is displayed in blue, with the setup employed in this paper (except for not including the hexadecapole), with the local Lagrangian bias relations of eq. (A.2) and fitted scales of $0.02 < k [hMpc^{-1}] < 0.15$. The horizontal and vertical dotted lines mark the expected values of the ShapeFit parameters. In the right panel of figure 3 we show the corresponding best-fit power spectrum multipoles (solid red and dotted blue lines) to the mean of the same 25 Abacus cubic mocks. The two models, which have been fitted with different binning and range of $k$, are almost indistinguishable, as expected. The figure highlights the excellent agreement between two-loop RPT and Folps with the minimal freedom setup, as well as their ability to recover unbiased values of the parameters. Both cases show minimal difference, well within the $1\sigma$ statistical error associated with a volume of 200 [Gpc h⁻¹]³.

We note that for the combination of parameters $m + n$ the Folps minimal freedom setup has a ~30% larger errorbars than the measurements from the RPT. This difference is due to the extra nuisance parameters that Folps considers (2 counter-terms and 1 extra stochastic term) that allow more freedom in the shape of the power spectrum. Nevertheless, the agreement between these two models in the best-fit value of the combined shape parameter is remarkable. The Folps maximal freedom setup presents larger errorbars, especially for the $m + n$ parameter due to the strong degeneracy between the $b_{3nl}$ and the shape parameters, as discussed already in [91].

These results indicate that both EFTofLSS and RPT models are able to reproduce with the same level of accuracy the ShapeFit cosmological constraints when the setup of nuisance parameters is set to be the same, reinforcing the robustness of the findings to different modeling choices: both types of approaches can be indistinctly employed to retrieve the cosmological constraints from DESI data. For this paper, we choose to use the RPT model as our baseline choice for describing the power spectrum multipoles.

### 3.3 Bispectrum modelling

Our bispectrum model has its foundations in the standard perturbation theory (SPT) formalism at tree-level, which has the form

$$B^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) = D_{\text{FoG}}^{B}(\mathbf{k}_1, \mathbf{k}_2) [2Z_1^{\text{SPT}}(\mathbf{k}_1)Z_1^{\text{SPT}}(\mathbf{k}_2)Z_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2)P_L(k_1)P_L(k_2) + 2\text{perm.}], \quad (3.6)$$

<footnote n="9">⁹Only for this plot we employ a slightly different version of the Abacus mocks, with a different HOD parametrization than those mocks described in section 2.1. The reason is that we aim to reproduce the results presented in the first batch of DESI full-shape comparison papers released in April 2024 [91, 96–98, 101] before the creation of the new set of Abacus mocks used in the rest of the paper.</footnote>

<footnote n="10">¹⁰This approach is slightly different from the one performed in the official DESI ShapeFit results, which is equivalent to setting the third-order bias $b_{3nl}$ to zero and leaving the tidal bias $b_{s2}$ free.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left panel: posteriors of the ShapeFit parameters obtained by fitting the power spectrum monopole and quadrupole of the 25 Abacus mocks for the LRG sample at $z = 0.8$, with a covariance corresponding to a volume of 200 [Gpc$^{-1}$]$^3$. The light orange and red contours display the performance of the EFTofLSS code Folps [91] in the range $0.02 < k [hMpc^{-1}] < 0.18$, for two types of galaxy bias setup: minimal freedom in filled contours and maximal freedom in empty contours, as labelled. The blue posteriors display the performance for the RPT model, which is used in this work for describing the power spectrum signal, and that has encoded a similar setup as the ‘minimal freedom’ for Folps, but with $0.02 < k [hMpc^{-1}] < 0.15$. The horizontal and vertical black dotted lines mark the expected values for each parameter given the true cosmology. Right panel: The upper and bottom sub-panels display the power spectrum monopole and quadrupole, respectively, for the Folps minimal freedom and RPT best-fitting models to the measurements of the 25 Abacus cubic mocks (not displayed for clarity). The grey dashed region is the error associated with one Abacus mock, corresponding to a volume of $8 [Gpc^{-1}]^3$.</figcaption></figure>

where ‘perm’ refers to cyclic permutations of the $k_1, k_2, k_3$ vectors. All terms are detailed in appendix A.

The GEO-FPT model was developed in [56] to increase the range of scales of validity of the tree-level bispectrum model, and thus extract additional non-linear information. In GEO-FPT, the linear matter power spectrum, $P_L$, of equation (3.6) is promoted to the non-linear matter power spectrum, $P_{NL}$. Furthermore, the perturbation theory $Z_2^{SPT}$ kernel is modified as follows:

$$Z_2^{\text{GEO}} = Z_2^{\text{SPT}} \times \left[ f_1 + f_2 \frac{\cos(\theta_{\text{med}})}{\cos(\theta_{\text{max}})} + f_3 \frac{\cos(\theta_{\text{min}})}{\cos(\theta_{\text{max}})} + f_4 \frac{A}{A_{\text{norm}}} + f_5 \frac{A^2}{A_{\text{norm}}^2} \right] \quad (3.7)$$

where the coefficients $f_1, \dots, f_5$ were calibrated on N-body simulations. This specific parametrisation was proposed after observing, in figure 1 of [56], that the SPT bispectrum residuals had a clear quadratic dependence on the area of the triangle, and a mild dependence on the flatness of the triangle (which can be obtained as combinations between cosines of the
minimum, intermediate and maximum external angles of the triangle). Refs. [56, 74] show that this model, with the coefficients $f_1, \dots, f_5$ fixed to the values obtained in their initial fit, is robust against change in the underlying cosmology. Therefore, these coefficients are not free parameters in our analysis, which results in a model for the galaxy bispectrum of the form:

$$B_{\text{gal}}^{\text{GEO}}(\mathbf{k}_1, \mathbf{k}_2) = P_{\text{NL}}(\mathbf{k}_1)P_{\text{NL}}(\mathbf{k}_2)\mathcal{Z}(\mathbf{k}_1, \mathbf{k}_2) + 2\text{perm.}, \quad (3.8)$$

where $\mathcal{Z}(\mathbf{k}_i, \mathbf{k}_j) \equiv D_{\text{FoG}}^{\text{B}}(\mathbf{k}_i, \mathbf{k}_j)Z_1^{\text{SPT}}(\mathbf{k}_i)Z_1^{\text{SPT}}(\mathbf{k}_j)Z_2^{\text{GEO}}(\mathbf{k}_i, \mathbf{k}_j)$ and $D_{\text{FoG}}^{\text{B}}$ is the Fingers-of-God damping factor for the bispectrum. More detail on these terms can be found in appendix A.

The non-linear matter power spectrum $P_{\text{NL}}$ is rescaled through the linear power spectrum scaling of equation (3.3). This accounts for the dependence of the bispectrum on the combined shape parameter $m+n$.

As for the other cosmological parameters $\alpha_\parallel, \alpha_\perp, f$, they modify the bispectrum model signal through equations (A.5), (A.10), (A.12) in appendix A and equation 2.11 in [56], while the $\sigma_{s8}$ parameter enters as a rescaling of the different terms of the 2-loop power spectrum model, scaling as $\sigma_{s8}^2 \propto P_L$.

The adopted range of scales for the bispectrum in all redshift bins matches the range in which the bispectrum model has been validated: $0.02 < k^B [h\text{Mpc}^{-1}] < 0.12$. The binning is the same as for the power spectrum, $\Delta k = 0.01 h\text{Mpc}^{-1}$.

For modelling the convolution effect of the survey geometry, we make use of the approximation of [45, 47], where the bispectrum kernels are not affected by the window.

In this approximation, only the power spectrum part of the bispectrum model equation (3.8) (i.e. $P_{\text{NL}}$) is affected by the window while the other terms are left unchanged. We obtain the windowed non-linear power spectrum $P_{\text{NL}}^W$ as in [82, 99], and compute the convolved bispectrum monopole $B_0^W$ as

$$B_0^W(k_1, k_2, k_3) = \int_{-1}^{1} d\mu_1 \int_{0}^{2\pi} d\phi \mathcal{Z}(\mathbf{k}_1, \mathbf{k}_2) P_{\text{NL}}^W(\mathbf{k}_1) P_{\text{NL}}^W(\mathbf{k}_2) + 2\text{perm.} \quad (3.9)$$

In the above expression, $\mu, \theta, \phi$ are the standard choice of angles in the multipole expansion that we use [102], which are detailed in appendix A.

Equation (3.9) is only valid for the bispectrum monopole. A generalization to bispectrum multipoles with the effect of a window should not follow this form, and a more careful treatment is needed.<fnref n="11" /> Consequently, in this paper we only consider the bispectrum monopole. This has limitations, since [56] showed that, in this particular model, the inclusion of the bispectrum quadrupoles results in a gain in both precision and accuracy on most cosmological parameters. We leave the development of an appropriate window treatment for the bispectrum quadrupoles, and thus the analysis of DESI data with bispectrum monopole and quadrupoles, for future work [79].

<footnote n="11">¹¹In our work we employ the bispectrum multipole expansion defined in [102]. In the case of the alternative base expansion presented in [103], valuable progress has been made in [104, 105] towards modelling the bispectrum window convolution. However, since the multipole-expansion bases of [102] and [103] do not have a direct correspondence between them, the work presented in [104] for the base of [103] is not necessarily applicable in our framework.</footnote>
## 3.4 Likelihood and covariance estimation

We sample the parameters posteriors via Markov chain Monte Carlo (MCMC), using the BRASS<fnref n="12" /> [59, 81] code, which follows the Metropolis-Hastings algorithm [106, 107]. We set the Gelman-Rubin convergence criterion to be $R - 1 < 0.01$.<fnref n="13" /> As shown in table 2, we specify broad, uniform priors in all parameters except for the deviations from Poissonian shot-noise, which we set to a normal prior centred at Poisson value, $N(1, 0.3^2)$. As motivated by [108], we estimate the full covariance matrix, including all off-diagonal terms. For each redshift bin and data-vector, the covariance is obtained from the 1000 independent realisations of the EZmock samples. Each power spectrum multipole data-vector has 13 elements, while the bispectrum monopole has 150 elements. Therefore, our largest data-vector, $P_{024} + B_0$, has a size of 189 elements, which corresponds to less than 20% of the number of simulations used to estimate the covariance matrix, resulting in a Hartlap factor [109] of ~0.8.

As reported in [19, 61], the Fourier-space covariance estimated from the 1000 EZmock has some potential inaccuracies. Because of this, refs. [19, 61] suggest to rescale the overall amplitude of the mock-estimated covariance by the factors displayed in table 3 of [61]. These factors were necessary to account for the mismatch between the EZmock and the DR1 covariances seen in configuration space. In this work, we assume the same rescaling factor that has been studied and validated for the power spectrum works as well for the bispectrum analysis. As we will discuss in section 5.2, we find these rescaling values to be valid for our analysis as well, with a potential systematic error involved with a mis-estimation of the covariance matrix accounting for less than ~20% of the statistical DR1 errors.

We consider the log-likelihood, $log L$, to be Gaussian, and apply the Sellentin-Heavens correction [110] to account for the error in the covariance introduced by estimating it from $n = 1000$ realisations,

$$ \log L = -\frac{n}{2} \log \left( 1 + \frac{(\mathbf{D}_{\text{meas.}} - \mathbf{D}_{\text{model}}) \text{Cov}_{\text{mocks}}^{-1} (\mathbf{D}_{\text{meas.}} - \mathbf{D}_{\text{model}})^T}{n-1} \right), \quad (3.10) $$

where $\mathbf{D}$ is the data-vector in each case, obtained either from measurements (meas.) or from the modelling (model); $n$ is the number of sample realisations used for estimating the covariance matrix, Cov.

## 4 Systematic error budget

In this section, we present the potential main sources of systematic error we have identified, and quantify their effect on parameter inference, by performing various tests and analyses of the Abacus mocks.

<footnote n="12">¹²https://github.com/hectorgil/Brass.</footnote>

<footnote n="13">¹³In all the MCMC runs, we first obtain a proposal parameter covariance through iterative runs where we sequentially increase the step size. The last run, with a step size of 1.9, has a burn-in of 10,000 steps for each walker (chain), and a total of 50,000 steps per walker. At that point, the code checks for convergence according to the Gelman-Rubin criterion, and in all cases it was found to be converged. With these settings, we found convergence to be obtained irrespective of the number of walkers, and in particular for the LRG redshift bins we used 4 walkers, while for the QSO we used 18 walkers.</footnote>
<table><thead><tr><th>Cosmological parameters</th><th>Prior range</th></tr></thead><tbody><tr><td>α<sub>iso</sub></td><td>U(0.7, 1.3)</td></tr><tr><td>α<sub>AP</sub></td><td>U(0.7, 1.3)</td></tr><tr><td>f</td><td>U(0.0, 20)</td></tr><tr><td>σ<sub>s8</sub></td><td>U(0.0, 10)</td></tr><tr><td>m + n</td><td>U(-1.0, 1.0)</td></tr></tbody></table>

<table><caption>**Table 2.** Priors for both the cosmological and nuisance parameters. We note as $U(x_{min}, x_{max})$ an uniform distribution between $x_{min}$ and $x_{max}$, and as $N(μ, σ^2)$ as a Gaussian distribution with mean $μ$ and standard deviation $σ$. In all cases (except from the parameters $A_P$, $A_B$, which have a Gaussian prior centred at 1 with $σ = 0.3$), the MCMC sampling is not affected by the boundaries of the priors on the remaining parameters thus being effectively improper priors.</caption><thead><tr><th>Nuisance parameters</th><th>Priors</th></tr></thead><tbody><tr><td>b<sub>1</sub></td><td>U(0, 10)</td></tr><tr><td>b<sub>2</sub></td><td>U(-20, 20)</td></tr><tr><td>A<sub>P</sub></td><td>N(1, 0.3)</td></tr><tr><td>A<sub>B</sub></td><td>N(1, 0.3)</td></tr><tr><td>σ<sub>P</sub></td><td>U(0, 10)</td></tr><tr><td>σ<sub>B</sub></td><td>U(0, 10)</td></tr></tbody></table>

We categorise systematic errors as originating from the inaccuracies in the theoretical modelling (section 4.1), from the fibre assignment procedure (section 4.2), from the impact of the HOD (section 4.3), and from the choice of fiducial cosmology (section 4.4).<fnref n="14" />

As in [61], we take into account a systematic error for a parameter once it is large enough to be correctly quantified on the mocks, which for DR1 translates into when the systematic shift is 20% or more of the DR1 statistical error. Following [61] the systematics are combined quadratically into a final (statistical plus systematic) error bar.

## 4.1 Theoretical systematics

For each redshift bin and galaxy sample, we consider the mean of the power spectrum and bispectrum measurements for all the 25 corresponding ‘complete’ Abacus mocks (see section 2.1), corresponding to the signal of a volume of 25 times the volume of DESI DR1. Consequently, we rescale the covariance matrix (obtained with the corresponding set of EZmocks) by a factor of 25. For all the data-vectors of interest in this work, {$P_{02}$, $P_{024}$, $P_{02} + B_0$, $P_{024} + B_0$}, we run the pipeline presented in section 3 on the mocks in the same way as we would do on the data.

Since they do not include any effects from observational sources of error, the ‘complete’ mocks can be seen as idealised simulations and we consider them as reference to assess the theoretical model systematics. These include both possible limitations or inaccuracies of

<footnote n="14">¹⁴The systematic error associated with the covariance matrix estimation is addressed by rescaling the covariance matrices by the factors given by table 3 of [61], which we check to be a valid approach in section 5.2.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Recovered cosmological parameters {$df, d\sigma_{s8}, \alpha_{\text{iso}}, \alpha_{\text{AP}}, m + n$} posteriors (68 and 95% joint confidence levels C.L.) for the joint power spectrum-bispectrum data-vector $P_{024} + B_0 = \{P_0, P_2, P_4, B_0\}$ from the mean of 25 ABACUS SUMMIT LRG complete mocks, and $P_{02} + B_0 = \{P_0, P_2, B_0\}$ from the mean of 25 ABACUS SUMMIT QSO complete mocks, which correspond to our baseline choice of multipoles. The panels of the figure display the posteriors for the four redshift bins considered. The employed covariance has been normalised by 25 in order to match the volume of the sample. Thus, the displayed posteriors correspond to ~25 times the DESI DR1 volume for each of the samples. We quantify the theoretical systematic error as the discrepancy between the maximum a posteriori values (MAP) and the truth for each parameter and galaxy sample.</figcaption></figure>
the power spectrum and bispectrum models, as well as effects of the analysis choices such as scale cuts, and the approximation introduced by the window convolution procedure. For each sample and cosmological parameter, the effect of the theoretical model systematics is given by the absolute difference between the true and the recovered values.<fnref n="15" />

In figure 4 we show the recovered parameter constraints from the complete Abacus mocks in the four redshift bins of study, compared with the true expected values (black dotted lines), for our baseline data-vector choices, $P_{024} + B_0$ for the LRG sample, and $P_{02} + B_0$ for the QSO sample. Instead of $f$ and $\sigma_{s8}$, to facilitate visualization, we display the posteriors in terms of $df \equiv f/f^{\text{true}}$, $d\sigma_{s8} \equiv \sigma_{s8}/\sigma_{s8}^{\text{true}}$. We quantify the residuals between the MAP values of each posterior and the underlying true values for the four parameters of interest in the rows labelled as 'Modelling' in tables 4, 5, 9, 10, for all data-vector combinations (not all of them are displayed in figure 4).

All tracers are consistent with the true $\alpha_{\text{iso}}$, $\alpha_{\text{AP}}$ at $2\sigma$ (corresponding to a volume ~25×DR1). The LRG2 redshift bin features a systematic shift in $d\sigma_{s8}$, while all LRGs have an offset in the $m+n$ parameter with respect to its expected value. We have previously seen that systematic shifts in $df$, $d\sigma_{s8}$ are anticipated when not including the bispectrum quadrupoles in the analysis [56]. As for $m+n$, we see in the systematic tables, which we refer to above, that including $B$ causes a slightly larger systematic shift in that parameter with respect to the power spectrum-only case. Given that $m+n$ is sensitive to the freedom of the bias model, as observed in figure 3, we leave for further work the investigation of the effect of assuming the local Lagrangian bias expansion (section 3.1) in a joint power spectrum-bispectrum analysis. Another possible source for this could be the approximation that we use for the bispectrum window convolution (see section 3), which is part of our current work.

## 4.2 Fibre assignment and imaging systematics

The DESI survey has two observational effects that can impact the clustering analysis: the fibre incompleteness and the photometric angular density fluctuations. Fibre incompleteness is caused by the inability of DESI to observe close pairs of objects within the patrol radius of a positioner. This is mitigated by performing several passes over the same area during the whole period of five years of observations. For this reason, this effect will be less important at the end of the five-year programme than for DR1. The angular density fluctuations are already present in the imaging maps, and are mitigated by the use of imaging weights in the measurement of the power spectra and bispectra, denoted imaging systematics. A thorough description of these can be found in refs. [19, 67, 111].

As for the imaging systematics, refs. [61, 67] found that in all LRG redshift bins the default linear imaging weights prescription is sufficient to capture the effect of observational conditions, whereas the QSO sample only suffers from a mild dependence on the systematic weights. The sample most affected by this effect is the ELG sample, which is the main reason not to include it in this work. For the QSO sample, we show in figure 6 the discrepancy

<footnote n="15">¹⁵We consider the *maximum a posteriori* (MAP) from the MCMC output as the recovered parameter value instead of the central value of the posterior distribution. The MAP represents the set of cosmological and nuisance parameters that maximise the posterior distribution. However, in most cases in our analysis (especially when fitting the mean of 25 mocks), the MAP and central value are very close.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of the MAP values for the cosmological parameters inferred from the complete and AltMTL mocks for all the tracers we consider in this work and for the choice of data-vectors of our baseline analysis. The dashed orange area is the $1-2\sigma$ region of the DESI DR1 analysis, that we have obtained from the blinded measurements in appendix B, while the individual error bars correspond to the parameter inference with the mean of the 25 mocks (thus rescaling the associated errors by 25). For the QSO sample we only report the product $df\sigma_{s8}$ as for this sample this degeneracy is not lifted (see section 3).</figcaption></figure>
(averaged over the 25 QSO AltMTL Abacus mocks) between the measurements with imaging weights included (e.g. $P_{\ell,weight}$) and without them (e.g. $P_{\ell,noweight}$), for all power spectrum multipoles and the bispectrum monopole. We observe the same trends as in figure 9 of [61] (which displays the results for ELG tracers) and that these results are consistent with the QSO imaging shifts quantified in [61, 67]. In [61] the correction is applied by modifying the theoretical power spectrum model $P_{\ell}^{\text{theory}}(k)$ to become an effective power spectrum $P_{\ell}^{\text{eff}}(k)$:

$$P_{\ell}^{\text{eff}}(k) = P_{\ell}^{\text{theory}}(k) + s_p \text{Poly}[(P_{\ell,weight}^{\text{data}} - P_{\ell,noweight}^{\text{data}}) - (P_{\ell,weight}^{\text{mock}} - P_{\ell,noweight}^{\text{mock}})] \\ + \text{Poly}[P_{\ell,weight}^{\text{mock}} - P_{\ell,noweight}^{\text{mock}}], \qquad (4.1)$$

and similarly for the bispectrum, where ‘Poly’ refers to a 3rd degree polynomial fit of the data points.<fnref n="16" /> The first correction term (varied with a nuisance parameter $s_p$) marginalises over the difference between imaging weights in mocks and in the DR1 data. The last term is a constant angular mode removal term.

Ref. [67] shows that the QSO sample is only mildly affected by the mode removal term, and marginalizing over the first term has no significant effect on the recovered constraints. The shifts presented in figure 6, denoted as ‘mode removal’ correction, correspond to the last term of equation (4.1). Therefore we do not include the marginalisation over the non-linearities of the imaging weights (the term modulated by $s_p$ in equation (4.1)). We leave for future work an implementation of this marginalization, which will enable us to include the bispectrum clustering of the ELG sample.

In figure 7 we show the effect of the mode removal correction in the blinded DR1 QSO data, by comparing the recovered constraints with and without the correction computed with the last term of equation (4.1). The effect of this correction is extremely small, both for the case shown in the figure ($P_{02} + B_0$) as for the other data-vector combinations. In all cases, the shifts in the posteriors are less than 20% of the corresponding DR1 errors, with part of the shift being absorbed by the $b_2$ and $\sigma_B$ parameters. In our analysis, we will model the QSO theoretical model with the mode removal correction applied, and leave the LRG tracers (which reported a smaller imaging systematic weight dependence than QSO in [61, 67]) uncorrected.

As for the fibre assignment systematic, we quantify it as the discrepancy between the recovered maximum a posteriori (MAP) values from the complete and AltMTL sets of 25 Abacus mocks. We present the shifts between the complete and AltMTL MAPs in the rows labelled as ‘Fibre assignment’ in tables 4, 5, 9, 10, for each of the four data-vector combinations considered in this work.

In figure 5, for the baseline choice of data-vector combination as labelled in the x-axis, we display the MAP values for all cosmological parameters of interest and the four samples, together with their recovered error (corresponding to the volume of the 25 simulations), compared with the DESI DR1 $1-2\sigma$ error bars (obtained from the blinded data, which is shown in appendix B). We include the parameter combination $f\sigma_{88}$, for easier comparison to the QSO sample results (see section 3).

<footnote n="16">¹⁶This polynomial fit is used instead of the data points to reduce the noise associated with the data, although its effect is minor. In the case of the bispectrum, due to its complexity, we do not use any polynomial fits to the data.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Mode removal residuals for the QSO sample for all data-vectors considered in this work. In the top row and lower left panel we show respectively the shifts for the power spectrum monopole, quadrupole and hexadecapole, with errorbars corresponding to the volume spanned by 25 Abacus cubic mocks, $V = 200$ (Gpc/h)$^3$. In the lower right plot we show the mode removal correction for the bispectrum, normalized by the measurement error bars corresponding to 25 Abacus mocks. Note that the volume considered in these mocks is significantly larger than the actual DESI DR1 dataset, which makes the mode-removal effect to barely affect the ShapeFit parameters inferred from the data, as shown in figure 7.</figcaption></figure>

## 4.3 HOD systematics

The systematic errors related to the HOD originate from the limitation of the bias model considered here (two-free parameters $b_1$ and $b_2$ with local-Lagrangian values for $b_{s2}$ and $b_{3nl}$), to describe the matter-halo-galaxy connection.

We start by considering a set of Abacus galaxy mocks with variations of the DESI HOD baseline, for both LRG and QSO samples. These alternative-HOD mocks were generated by the DESI collaboration by considering variations about a baseline HOD chosen for the so-called second-generation mocks, which were calibrated with the DESI Early Data Release [16].

We perform our analysis on cubic box simulations spanning the baseline (standard) and the alternatives to the standard DESI HOD, which consist in variations in the modelling of velocity bias, environment-based bias and the masses of the halos hosting the central
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Cosmological parameters posteriors recovered from the blinded DR1 QSO catalogue, using the combination of data-vectors $P_{02} + B_0$, with (in purple) and without (in orange) the mode removal correction defined in equation (4.1). The effect of this correction is very mild, with the contours being effectively indistinguishable for DR1 statistical error bars.</figcaption></figure>

and satellite galaxies. We employ the same notation as in refs. [61, 112, 113] and refer the reader to section 2 in [112] for their description. As in these works, we consider seven alternatives (A1, A2, A3, B0, B1, B2, B3) to the standard (A0) HOD for LRG tracers, and 3 alternatives to the standard ($QSO_0$) HOD ($QSO_1, QSO_2, QSO_3$) for QSOs. We display the recovered posteriors of the LRG galaxies at $z = 0.8$ in figure 8, where the different colours correspond to different variations of HOD.<fnref n="17" />

<footnote n="17">¹⁷The comparison of constraints for the QSO alternative HODs is qualitatively similar to the LRGs displayed in figure 8, hence we do not explicitly show the analogous plot for QSOs.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Recovered posteriors (68 and 95 % joint C.L.) of cosmological parameters {$\alpha_{\text{iso}}, \alpha_{\text{AP}}, m + n, df, d\sigma_{\text{s8}}$} for the 8 alternative HOD variations (see [112] for a description), compared to the baseline HOD, noted as A0. The $P_{024} + B_0$ data-vector is obtained from each set of 25 Abacus HOD simulations, in cubic boxes describing LRG galaxies at $z = 0.8$, so the covariance is normalised by 25.</figcaption></figure>

We then estimate the contribution from the HOD systematic as the mean across HOD variations of the offset between their inferred MAP parameters and the MAP obtained from the standard HOD simulations (labelled as A0 in figure 8). We see in figure 8 that the $\alpha_{\text{iso}}, \alpha_{\text{AP}}, m + n$ parameters do not suffer significantly from HOD-related systematics in the LRG tracers, while there is a small error budget contribution for the $f$ and $\sigma_{\text{s8}}$ parameters. The values for these systematic offsets, for both the LRG tracers and the QSO are listed in the row labelled as ‘HOD’ in tables 4, 5, 9, 10, for each of the four data-vector combinations considered.
<table><caption>**Table 3.** Summary of the cosmologies considered to test the systematic error associated with the choice of fiducial cosmology (analogously as in [114, 116]). c000 is the baseline fiducial cosmology.</caption><thead><tr><th>Name</th><th>ω<sub>b</sub></th><th>ω<sub>cdm</sub></th><th>h</th><th>10<sup>9</sup>A<sub>s</sub></th><th>n<sub>s</sub></th><th>N<sub>ur</sub></th><th>w<sub>0</sub></th><th>w<sub>a</sub></th></tr></thead><tbody><tr><td>Planck ΛCDM – c000</td><td>0.02237</td><td>0.1200</td><td>0.6736</td><td>2.0830</td><td>0.9649</td><td>2.0328</td><td>-1</td><td>0</td></tr><tr><td>Low-Ω<sub>m</sub> – c001</td><td>0.02237</td><td>0.1134</td><td>0.7030</td><td>2.0376</td><td>0.9638</td><td>2.0328</td><td>-1</td><td>0</td></tr><tr><td>Thawing DE – c002</td><td>0.02237</td><td>0.1200</td><td>0.6278</td><td>2.3140</td><td>0.9649</td><td>2.0328</td><td>-0.7</td><td>-0.5</td></tr><tr><td>High-N<sub>eff</sub> – c003</td><td>0.02260</td><td>0.1291</td><td>0.7160</td><td>2.2438</td><td>0.9876</td><td>2.6868</td><td>-1</td><td>0</td></tr><tr><td>Low-σ<sub>8</sub> – c004</td><td>0.02237</td><td>0.1200</td><td>0.6736</td><td>1.7949</td><td>0.9649</td><td>2.0328</td><td>-1</td><td>0</td></tr></tbody></table>

## 4.4 Systematics due to fiducial cosmology assumption

The last source of systematic error that we explore is the effect of choice of fiducial cosmology. The fiducial cosmology that we use in our analysis, which is the one corresponding to the Abacus base simulations (c000, detailed in section 2.1, first line of table 3) is assumed in two different steps:

1) In the conversion from redshifts into distances in our galaxy and quasar catalogues.

2) In the calculation of our power spectrum linear matter template.

We keep fixed the power spectrum template, for which we pre-compute the perturbation theory integral terms through the PTCOOL<fnref n="18" /> code [45, 57], where the dependence on $\sigma_{s8}$ and $m+n$ is factorised out of the integrals. This approximation, which was validated in appendix D of [60], saves significant time and resources by eliminating the need to re-compute the perturbation theory integrals in each step of the MCMC.

Consequently, we need to quantify the effect of assuming a fiducial cosmology that could be different from the true cosmology. In this work, we only focus on the systematic error coming from the choice of fiducial cosmology in the template (i.e., in point 2). We neglect the impact of the fiducial cosmology in the catalogue creation, which has been estimated to be negligible in [114, 115].

Following the approach taken by the DESI collaboration in the main BAO and Full-Shape ShapeFit analyses [114, 116], we consider four alternative fiducial cosmologies, which are summarised in table 3. They correspond to the first four alternative cosmologies of the Abacus suite of simulations, which are: c001 (low-$\Omega_m$), c002 (thawing dark energy), c003 (high-$N_{\text{eff}}$), c004 (low-$\sigma_8$).

The dilation parameters, $\alpha_\parallel$, $\alpha_\perp$, and the shape parameter $m+n$, are sensitive to mismatches between fiducial and true cosmology. This means that when the fiducial cosmology is different from the true one, we expect $\alpha_\parallel, \alpha_\perp$ to differ from 1, and $m+n$ to be non-zero, but knowing the two cosmologies we can compute their expected values.

For a given parameter $x$, we define the variable $\Delta x = x - x^{\text{exp}}$ as the shift of the obtained parameter with respect to its expected value given the true and fiducial cosmologies. For $f, \sigma_{s8}$ and $f\sigma_{s8}$ we report their obtained value, since these parameters are not expected to depend on the choice of template. Therefore, any deviation in the recovered value from

<footnote n="18">¹⁸https://github.com/hectorgil/PTcool.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Similar as figure 8 but showing the differences between the parameter constraints (for $P_{024} + B_0$ of Abacus LRGs at $z=0.8$) while assuming a template based on each of the alternative fiducial cosmologies: c001, c002, c003, c004. Given that these fiducial cosmologies differ from the true one of the simulations, $\Delta$ indicates the discrepancy between the recovered and expected parameters. The posteriors have been obtained from the mean of the 25 LRG Abacus cubic boxes at $z = 0.8$.</figcaption></figure>

the expected parameter value when assuming a different fiducial cosmology is a signal of a systematic error and must be quantified.

In figure 9 we display the residuals between the recovered cosmological parameters and their expected values, when using the power spectrum templates for each of the alternative fiducial cosmologies, for LRG tracers in the average of the 25 respective cubic boxes at $z = 0.8$. If there was no systematic error induced by the choice of fiducial cosmology, all the posterior distributions in figure 9 would be coincident. The effect of the fiducial cosmology systematics
<table><caption>**Table 4.** Contributions to the systematic error budget for each parameter and redshift bin, for the combination of summary statistics of $P_{024} + B_0$, i.e. the power spectrum monopole, quadrupole and hexadecapole, plus the bispectrum monopole. All errors are expressed here as percentage of the DR1 $1\sigma$ statistical error. The threshold to propagate a specific contribution to the systematic error budget is 20% of the corresponding DR1 $1\sigma$ statistical error (see text for more details). For the HOD and fiducial cosmology systematics, we also show in parenthesis the maximum offsets obtained among all HOD simulations and among all fiducial cosmology choices, respectively. The last row shows the quadratic sum of all the contributing systematic sources, which is to be added in quadrature to the statistical error bars obtained from the corresponding cosmological parameter posteriors.</caption><thead><tr><th>[%σ<sub>DR1</sub>]</th><th>P<sub>024</sub> + B<sub>0</sub></th><th>σ<sub>α<sub>iso</sub></sub></th><th>σ<sub>α<sub>AP</sub></sub></th><th>σ<sub>f</sub>σ<sub>s8</sub></th><th>σ<sub>m+n</sub></th><th>σ<sub>f</sub></th><th>σ<sub>σ<sub>s8</sub></sub></th></tr></thead><tbody><tr><td rowspan="4">Modelling</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 20</td><td>~ 50</td><td>&lt; 20</td><td>&lt; 10</td></tr><tr><td>LRG2</td><td>&lt; 10</td><td>~ 20</td><td>~ 50</td><td>~ 50</td><td>~ 40</td><td>~ 90</td></tr><tr><td>LRG3</td><td>~ 20</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 50</td><td>~ 40</td><td>~ 40</td></tr><tr><td>QSO</td><td>~ 20</td><td>&lt; 10</td><td>~ 40</td><td>&lt; 20</td><td>—</td><td>—</td></tr><tr><td rowspan="4">Observational</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 20</td><td>~ 30</td></tr><tr><td>LRG2</td><td>&lt; 10</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 30</td><td>~ 50</td><td>~ 70</td></tr><tr><td>LRG3</td><td>&lt; 20</td><td>~ 50</td><td>~ 60</td><td>&lt; 10</td><td>&lt; 10</td><td>~ 50</td></tr><tr><td>QSO</td><td>~ 40</td><td>~ 40</td><td>~ 90</td><td>&lt; 20</td><td>—</td><td>—</td></tr><tr><td rowspan="4">HOD</td><td>LRG1</td><td>&lt; 10(16)</td><td>&lt; 10(14)</td><td>&lt; 10(8)</td><td>&lt; 10(8)</td><td>&lt; 10(30)</td><td>&lt; 20(39)</td></tr><tr><td>LRG2</td><td>&lt; 10(16)</td><td>&lt; 10(13)</td><td>&lt; 10(10)</td><td>&lt; 10(9)</td><td>&lt; 10(24)</td><td>&lt; 20(33)</td></tr><tr><td>LRG3</td><td>&lt; 10(20)</td><td>&lt; 10(16)</td><td>&lt; 10(10)</td><td>&lt; 10(8)</td><td>&lt; 10(19)</td><td>&lt; 20(42)</td></tr><tr><td>QSO</td><td>&lt; 10(19)</td><td>&lt; 10(22)</td><td>&lt; 20(39)</td><td>&lt; 20(39)</td><td>—</td><td>—</td></tr><tr><td rowspan="4">Fiducial cosmology</td><td>LRG1</td><td>&lt; 10(15)</td><td>&lt; 10(10)</td><td>&lt; 20(15)</td><td>&lt; 20(70)</td><td>~ 30(62)</td><td>~ 30(69)</td></tr><tr><td>LRG2</td><td>&lt; 10(15)</td><td>&lt; 10(10)</td><td>&lt; 20(18)</td><td>&lt; 20(77)</td><td>~ 30(50)</td><td>~ 20(59)</td></tr><tr><td>LRG3</td><td>&lt; 10(19)</td><td>&lt; 10(12)</td><td>&lt; 20(19)</td><td>&lt; 20(75)</td><td>~ 20(39)</td><td>~ 30(75)</td></tr><tr><td>QSO</td><td>&lt; 10(23)</td><td>&lt; 10(18)</td><td>&lt; 10(19)</td><td>~ 30(71)</td><td>—</td><td>—</td></tr><tr><td rowspan="4">Total</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 20</td><td>~ 50</td><td>~ 30</td><td>~ 42</td></tr><tr><td>LRG2</td><td>&lt; 20</td><td>~ 20</td><td>~ 50</td><td>~ 58</td><td>~ 71</td><td>~ 116</td></tr><tr><td>LRG3</td><td>~ 20</td><td>~ 50</td><td>~ 60</td><td>~ 50</td><td>~ 45</td><td>~ 71</td></tr><tr><td>QSO</td><td>~ 45</td><td>~ 40</td><td>~ 103</td><td>~ 30</td><td>—</td><td>—</td></tr></tbody></table>

is quantified as the mean across samples of the shifts between these residuals. As seen in the row labelled as 'Fiducial' in tables 4, 5, 9 and 10, the contribution of the choice of fiducial cosmology to the systematic error budget is generally minor: the main sources of systematic error come from the fibre assignment and modelling parts of the analysis.

## 5 Results

While this analysis is not required by the DESI collaboration rules to adhere to the DESI blinding procedure (the data has been unblinded for several months at the time of writing this paper), we opt to work with blinded data (as explained in section 3) until we freeze the analysis choices and obtain our systematic error budget, similarly as is done in the DESI
<table><caption>**Table 5.** Same as table 4, but for the power spectrum monopole and quadrupole, plus the bispectrum monopole ($P_{02} + B_0$).</caption><thead><tr><th>[%σDR1]</th><th>P<sub>02</sub> + B<sub>0</sub></th><th>σ<sub>αiso</sub></th><th>σ<sub>αAP</sub></th><th>σ<sub>f</sub>σ<sub>s8</sub></th><th>σ<sub>m+n</sub></th><th>σ<sub>f</sub></th><th>σ<sub>σs8</sub></th></tr></thead><tbody><tr><td rowspan="4">Modelling</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 20</td><td>~ 50</td><td>&lt; 20</td><td>&lt; 20</td></tr><tr><td>LRG2</td><td>&lt; 10</td><td>&lt; 10</td><td>~ 30</td><td>~ 50</td><td>~ 70</td><td>~ 70</td></tr><tr><td>LRG3</td><td>~ 20</td><td>&lt; 20</td><td>~ 30</td><td>~ 50</td><td>~ 90</td><td>~ 70</td></tr><tr><td>QSO</td><td>~ 30</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 20</td><td>—</td><td>—</td></tr><tr><td rowspan="4">Observational</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 20</td><td>~ 40</td></tr><tr><td>LRG2</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 10</td><td>~ 40</td><td>~ 60</td><td>~ 50</td></tr><tr><td>LRG3</td><td>&lt; 20</td><td>&lt; 10</td><td>~ 30</td><td>&lt; 10</td><td>~ 60</td><td>~ 90</td></tr><tr><td>QSO</td><td>~ 40</td><td>&lt; 20</td><td>~ 50</td><td>&lt; 20</td><td>—</td><td>—</td></tr><tr><td rowspan="4">HOD</td><td>LRG1</td><td>&lt; 10(16)</td><td>&lt; 10(12)</td><td>&lt; 10(7)</td><td>&lt; 10(11)</td><td>&lt; 10(27)</td><td>&lt; 10(36)</td></tr><tr><td>LRG2</td><td>&lt; 10(15)</td><td>&lt; 10(8)</td><td>&lt; 10(7)</td><td>&lt; 10(11)</td><td>&lt; 20(21)</td><td>&lt; 10(25)</td></tr><tr><td>LRG3</td><td>&lt; 10(19)</td><td>&lt; 10(11)</td><td>&lt; 10(8)</td><td>&lt; 10(11)</td><td>&lt; 10(21)</td><td>&lt; 20(36)</td></tr><tr><td>QSO</td><td>&lt; 10(12)</td><td>&lt; 10(7)</td><td>&lt; 20(14)</td><td>&lt; 20(26)</td><td>—</td><td>—</td></tr><tr><td rowspan="4">Fiducial cosmology</td><td>LRG1</td><td>&lt; 10(19)</td><td>&lt; 10(18)</td><td>&lt; 20(26)</td><td>&lt; 20(64)</td><td>~ 30(70)</td><td>~ 30(70)</td></tr><tr><td>LRG2</td><td>&lt; 10(19)</td><td>&lt; 10(17)</td><td>&lt; 20(30)</td><td>&lt; 20(70)</td><td>~ 30(57)</td><td>~ 20(59)</td></tr><tr><td>LRG3</td><td>&lt; 10(23)</td><td>&lt; 10(21)</td><td>&lt; 20(33)</td><td>&lt; 20(68)</td><td>~ 30(44)</td><td>~ 30(75)</td></tr><tr><td>QSO</td><td>&lt; 20(36)</td><td>&lt; 10(14)</td><td>&lt; 10(18)</td><td>~ 20(66)</td><td>—</td><td>—</td></tr><tr><td rowspan="4">Total</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 20</td><td>~ 50</td><td>~ 30</td><td>~ 50</td></tr><tr><td>LRG2</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 30</td><td>~ 64</td><td>~ 97</td><td>~ 88</td></tr><tr><td>LRG3</td><td>~ 20</td><td>&lt; 20</td><td>~ 42</td><td>~ 50</td><td>~ 112</td><td>~ 118</td></tr><tr><td>QSO</td><td>~ 50</td><td>&lt; 20</td><td>~ 50</td><td>~ 20</td><td>—</td><td>—</td></tr></tbody></table>

key papers [17, 61, 66]. We perform all of the tests that involve data on the official blinded data of DR1 (see section 2 and references therein). We report the results for the blinded data in appendix B, and for the true, unblinded data in section 5.1.

Considering the associated errors (as quantified on the mocks, section 4) and projection effects in blinded data (presented in appendix B) we choose our bispectrum baseline analysis to consist of the data-vector $P_{024} + B_0$ for the LRG1, LRG2, LRG3 redshift bins, and $P_{02} + B_0$ for the QSOs. In fact, tables 4, 5, 9, 10 and appendix B show that the inclusion of the power spectrum hexadecapole for the LRGs reduces the statistical errors without introducing additional systematics, but the inclusion of the power spectrum hexadecapole for the QSO increases notably the systematic error budget without reducing the statistical errors in any appreciable way. Hereafter this baseline combination is referred to as $P + B$.

The multipole combination for the power spectrum represents another difference with the official DESI analysis, [61], where the power spectrum hexadecapole is not included. For comparison, we also present power spectrum only results, where the data-vectors involved are $P_{024}$ for the LRG1, LRG2, LRG3 redshift bins, and $P_{02}$ for the QSOs.
## 5.1 Unblinded data

After having frozen the pipeline, analyses choices and combinations of data-vectors, we perform the analysis on the measurements of the unblinded power spectra and bispectra baseline combination $P+B$.

Figure 10 presents the results for the baseline $P+B$ combination on the DESI DR1 data (see appendix B for detailed justification of the choice). To highlight the role of the bispectrum, the figure also displays the power spectrum-only results ($P_{024}$ for the LRG bins, and $P_{02}$ for the QSO, labelled as $P$), in addition to the official DESI results presented in ref. [61], in terms of ShapeFit parametrization labelled as 'DESI SF'. All constraints are highly compatible.

Despite the analysis differences (sample, selection, k-cuts, multipoles considered, different theoretical modelling), the $P$ results presented here are very consistent with the official DESI ones, highlighting the robustness of the analyses. Compared to $P$ only, the full $P+B$ analysis provides mild improvements on the size of the errors for the $\alpha_{\text{iso}}$ (9% improvement), $f\sigma_{\text{s8}}$ (9% improvement) and $m+n$ (11% improvement) parameters. As expected, since the bispectrum monopole has no significant anisotropic signal, there is no change in the size of the constraints for the parameter $\alpha_{\text{AP}}$. The detailed values of the maximum a posteriori (MAP) and the total error budget (including statistical and systematic contributions) are reported in table 6. The corresponding nuisance parameters are reported in appendix C.

Notably, the inclusion of the bispectrum monopole breaks the $f-\sigma_{\text{s8}}$ degeneracy. We denote with the prefix ‘d’ the ratio of $f$ or $\sigma_{\text{s8}}$ with respect to its fiducial (c000 cosmology) value at any effective redshift. The fiducial values are: $f^{\text{fid}} = \{0.763, 0.817, 0.870, 0.928\}$ and $\sigma_{\text{s8}}^{\text{fid}} = \{0.621, 0.565, 0.501, 0.401\}$, corresponding to LRG1, LRG2, LRG3 and QSO. We obtain $df = \{0.888_{-0.089}^{+0.186}, 0.977_{-0.220}^{+0.182}, 1.030_{-0.085}^{+0.368}\}$, $d\sigma_{\text{s8}} = \{1.224_{-0.133}^{+0.091}, 1.071_{-0.163}^{+0.278}, 1.000_{-0.223}^{+0.088}\}$ respectively for the LRG1, LRG2, LRG3 bins. This is the main result of this work.

The combination of all the redshift bins yield a cumulative 10.1% constraint on $f$ (with less than 10% of the contribution to the total error bar coming from the systematic errors), and of 8.4% on $\sigma_{\text{s8}}$ (27% of the total error bars is the contribution from systematics). For the other compressed parameters the cumulative errors are $\sigma_{\alpha_{\text{iso}}}$ = 0.9% (9% improvement with respect to our power spectrum-only analysis); $\sigma_{\alpha_{\text{AP}}}$ = 2.3% (no improvement with respect to power spectrum-only analysis); $\sigma_{f\sigma_{\text{s8}}}$ = 5.1% (9% improvement); $\sigma_{m+n}$ = 2.3% (11% improvement). For these parameters, the contribution to the cumulative error from the systematics is always below 11% of the total — for the $\alpha$ parameters is below 3%. A discussion on the $\chi^2$ values can be found in section 5.2.

Figure 11 provides an alternative visualization of the $P+B$ results (solid contours) as well as the $P$-only results (dashed contours), expressed in terms of the evolution across the samples (in different colours as labelled) of the angular diameter distance $D_{\text{M}}/r_{\text{d}}$ and Hubble distance $D_{\text{H}}/r_{\text{d}}$, together with the $f\sigma_{\text{s8}}$ and $m+n$ combined shape parameter.

It is interesting to compare figure 11 with figure 2 of Reference [81], which displays the analogous constraints for BOSS and eBOSS when only employing the power spectrum multipoles of the LRG and QSO samples; and also including the BAO-only Ly$\alpha$ results (which we do not consider in this work). It is remarkable that, with only one year of observations, DESI data yield comparable constraints. Subsequent DESI data releases are expected to provide significant improvements.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Constraints (68% C.L.) on the parameters {$\alpha_{\text{iso}}$, $\alpha_{\text{AP}}$, $df\sigma_{\text{s8}}$, $m+n$, $df$, $d\sigma_{\text{s8}}$} for the unblinded DESI data in the four redshift bins (LRG1, LRG2, LRG3, QSO) as indicated in abscissa. The different colours correspond to the baseline $P+B$ (orange), $P$ only (purple) and the official DESI SF (ShapeFit) analysis [61] results (blue). The dashed line marks the fiducial, c000, cosmology, and $df$, $d\sigma_{\text{s8}}$, and $df\sigma_{\text{s8}}$ are computed with respect to the fiducial model (i.e., $f/f^{\text{c000}}$, $\sigma_{\text{s8}}/\sigma_{\text{s8}}^{\text{c000}}$). Despite the analysis differences, the $P$ results presented here are very consistent with the official DESI ones. The addition of the bispectrum breaks the $f\sigma_{\text{s8}}$ degeneracy, in addition it tightens the error bars especially on the $df\sigma_{\text{s8}}$, $\alpha_{\text{iso}}$ and $m+n$ parameters. A rigorous interpretation of the ShapeFit constraints in terms of ΛCDM and extensions of it will be provided in [79].</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Posterior distributions (68 and 95% C.L.) for the parameters {$D_H(z_{\text{eff}})/r_d, D_M(z_{\text{eff}})/r_d, f\sigma_{s8}(z_{\text{eff}}), m(z_{\text{eff}}) + n(z_{\text{eff}})$} in the four redshift bins for the baseline analysis ($P+B$). The theoretical evolution of the cosmological parameters according to Planck 2018 ΛCDM (corresponding to the fiducial c000 cosmology) is shown as a black line. The 2D dashed contours represent the power spectrum-only analysis highlighting how the inclusion of the bispectrum improves significantly the constraints on the $m(z) + n(z)$ combined shape parameter. The effective redshift, $z_{\text{eff}}$, is defined in equation (2.1).</figcaption></figure>
<table><caption>**Table 6.** Results from our two main analyses, $P+B$ (including the bispectrum monopole) and $P$ (only power spectrum), for the LRG and QSO redshift bins. In each case, we show the maximum a posteriori (MAP), together with the $1\sigma$ region (which accounts for both statistical and systematic errors) centred in the MAP value, for the ShapeFit parameters $\alpha_{iso}$, $\alpha_{AP}$, $f\sigma_{s8}$, $m+n$, together with the parameters $f$ and $\sigma_{s8}$. The prefix ‘d’ in a parameter indicates its ratio with respect to its fiducial value at the effective redshift of the bin. The corresponding fiducial values are: $f^{\text{fid}} = \{0.763, 0.817, 0.870, 0.928\}$ and $\sigma_{s8}^{\text{fid}} = \{0.621, 0.565, 0.501, 0.401\}$, for respectively LRG1, LRG2, LRG3 and QSO. We additionally report the corresponding values for $\chi^2$ over the degrees of freedom for each case (see section 5.2 for discussion).</caption><thead><tr><th></th><th>Sample</th><th>α<sub>iso</sub></th><th>α<sub>AP</sub></th><th>dfσ<sub>s8</sub></th><th>m + n</th><th>df</th><th>dσ<sub>s8</sub></th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td rowspan="4">P + B</td><td>LRG1</td><td>0.989<sup>+0.031</sup><sub>-0.007</sub></td><td>1.025<sup>+0.044</sup><sub>-0.043</sub></td><td>1.087<sup>+0.133</sup><sub>-0.075</sub></td><td>-0.025<sup>+0.044</sup><sub>-0.053</sub></td><td>0.888<sup>+0.186</sup><sub>-0.089</sub></td><td>1.224<sup>+0.091</sup><sub>-0.133</sub></td><td>241/178</td></tr><tr><td>LRG2</td><td>0.962<sup>+0.025</sup><sub>-0.012</sub></td><td>1.012<sup>+0.052</sup><sub>-0.038</sub></td><td>1.047<sup>+0.099</sup><sub>-0.095</sub></td><td>0.042<sup>+0.023</sup><sub>-0.070</sub></td><td>0.977<sup>+0.182</sup><sub>-0.220</sub></td><td>1.071<sup>+0.278</sup><sub>-0.163</sub></td><td>291/178</td></tr><tr><td>LRG3</td><td>0.990<sup>+0.025</sup><sub>-0.006</sub></td><td>0.967<sup>+0.046</sup><sub>-0.036</sub></td><td>1.029<sup>+0.137</sup><sub>-0.062</sub></td><td>0.032<sup>+0.024</sup><sub>-0.068</sub></td><td>1.030<sup>+0.368</sup><sub>-0.085</sub></td><td>1.000<sup>+0.088</sup><sub>-0.223</sub></td><td>253/178</td></tr><tr><td>QSO</td><td>1.020<sup>+0.021</sup><sub>-0.025</sub></td><td>1.000<sup>+0.082</sup><sub>-0.053</sub></td><td>1.113<sup>+0.086</sup><sub>-0.126</sub></td><td>0.065<sup>+0.035</sup><sub>-0.046</sub></td><td>—</td><td>—</td><td>193/166</td></tr><tr><td rowspan="4">P</td><td>LRG1</td><td>0.964<sup>+0.027</sup><sub>-0.013</sub></td><td>1.024<sup>+0.059</sup><sub>-0.030</sub></td><td>1.001<sup>+0.136</sup><sub>-0.066</sub></td><td>0.003<sup>+0.020</sup><sub>-0.084</sub></td><td>—</td><td>—</td><td>33/31</td></tr><tr><td>LRG2</td><td>0.955<sup>+0.012</sup><sub>-0.023</sub></td><td>1.012<sup>+0.049</sup><sub>-0.035</sub></td><td>1.010<sup>+0.128</sup><sub>-0.073</sub></td><td>0.041<sup>+0.016</sup><sub>-0.093</sub></td><td>—</td><td>—</td><td>37/31</td></tr><tr><td>LRG3</td><td>0.988<sup>+0.020</sup><sub>-0.019</sub></td><td>0.935<sup>+0.045</sup><sub>-0.035</sub></td><td>1.111<sup>+0.146</sup><sub>-0.095</sub></td><td>-0.004<sup>+0.040</sup><sub>-0.073</sub></td><td>—</td><td>—</td><td>31/31</td></tr><tr><td>QSO</td><td>1.006<sup>+0.031</sup><sub>-0.022</sub></td><td>1.010<sup>+0.081</sup><sub>-0.061</sub></td><td>1.105<sup>+0.174</sup><sub>-0.088</sub></td><td>0.067<sup>+0.021</sup><sub>-0.069</sub></td><td>—</td><td>—</td><td>20/18</td></tr></tbody></table>

Within the considered redshift bins, our baseline $P+B$ analysis shows an average improvement with respect to the DESI SF results [61] of ~10% for $\alpha_{AP}$; of ~20% for $f\sigma_{s8}$; and of ~30% for $m+n$. We obtain ~5% larger constraints than DESI SF for the parameter $\alpha_{iso}$. Analogously, our $P$-only analysis results in an average improvement in the coincident redshift bins (LRG1, LRG2, LRG3, QSO) with respect to DESI SF [61] of ~9% for $\alpha_{AP}$; of ~4% for $f\sigma_{s8}$; of ~15% for $m+n$; and ~15% larger constraints than DESI SF for $\alpha_{iso}$.<fnref n="19" />

The interpretation of the ShapeFit compressed parameters presented in this work in terms of constraints on cosmological parameters of a specific model, such as the primordial size of fluctuations, $A_s$, the matter density, $\Omega_m$, and the expansion parameter $H_0$; as well as the comparison with CMB data will be reported in ref. [79].

## 5.2 Statistical error validation

We aim to validate the statistical errors component obtained on the unblinded data with the results obtained from Abacus AltMTL mocks. The reasons are two-fold.

1. To check whether the DESI DR1 data catalogues have similar statistical properties to the Abacus mocks used to validate the modelling and to test for potential systematic errors.

2. To understand the reason for the high minimum-$\chi^2$ values found when fitting the power spectrum and bispectrum from the LRG data (reported in table 6); and to quantify the potential impact of this high minimum-$\chi^2$ values on the determination of the statistical error component.

<footnote n="19">¹⁹The fact that the $\alpha_{iso}$ parameter is slightly less constrained in our power spectrum approach is probably due to having a lower $k_{\text{max}}$.</footnote>
3. To validate that the rescaling of the covariance that we have applied, based on section 5.7 of ref. [61] is still applicable within our analysis involving the power spectrum hexadecapole and the bispectrum monopole signals.

In order to do so, we perform the fit on the 25 individual realisations of the LRG and QSO Abacus AltMTL mocks, which contain the DESI DR1 data survey geometry and mimic the fibre assignment effect. We employ the same covariance used for fitting the DR1 data, derived from the 1000 FFA EZmocks realisations, including the correction factors of table 3 of ref. [61] as explained in section 2.1. We do not include any extra systematic contribution, as here we are only interested in testing the relative statistical properties of the data compared to the mocks.

In figure 12 we display the distribution of the $1\sigma$ (or 68% CL) errors for the four parameters of interest, $\{\alpha_\parallel, \alpha_\perp, f\sigma_{s8}, m+n\}$, for the four redshift bins studied in this work: top-left panel for the LRG1; top-right panel for the LRG2; bottom-left panel for the LRG3; bottom-right panel for the QSO, as labelled. For each of these panels, the blue dots display the statistical errors of the 25 individual AltMTL Abacus mocks. The statistical errors corresponding to the DESI DR1 data are shown as the red crosses. We compare the distribution of statistical errors with the standard deviation computed from the medians of the 25 individual posteriors, shown in orange.<fnref n="20" /> In addition, the histogram of the best-fitting $\chi^2$ per number of degrees of freedom (dof) is also shown for each of the LRG and QSO samples, where the red vertical line represents the minimum $\chi^2$ found for the data (and reported in table 6), and the blue histogram is from the 25 mocks.

Figure 12 shows that the errors obtained from the data are in general very consistent with the typical errors we obtain from the 25 mocks. Only for the $m+n$ parameter, for the LRG3 and QSO bins, we find that the error of the data is at the high end of the distribution of errors from the mocks. We do not consider this outlier as something statistically worrisome, as we estimate that such outliers sit in a region of ~2$\sigma$ from the typical size of errors.

We also notice that for the LRG1 and the LRG3 bins, the distribution of $\chi^2$/dof of the Abacus AltMTL mocks is centred around values ~1.4 to 1.5 while for LRG2 and QSO samples it is centred closer to the expected value of 1.0. One possible explanation is that the model is not able to accurately describe the mocks. In some cases (see table 4 for the LRGs and table 5 for the QSO) we have reported modelling systematic errors of 20% to 50% of the total statistical error budget, that could explain the high $\chi^2$ values in some of these bins. Another possibility is that the mock-based covariance we are employing is slightly underestimating the true covariance of the mocks, even after the correction term of table 3 of ref. [61]. We leave for future work a more thorough exploration of the covariance of the bispectrum, as it requires the development of more realistic mocks and more sophisticated techniques for an analytical-based covariance, which goes beyond the scope of this paper.

Finally, we also highlight that the minimum $\chi^2$ found for the data is broadly in agreement with those obtained for the mocks, with maybe the only exception of the LRG2 bin, where the highest $\chi^2$/dof found on a mock is of ~1.4, whereas for the data is ~1.6. However, this

<footnote n="20">²⁰In the limit of a large number of realisations, the two estimates should coincide for a Gaussian distribution: within the limits of using only 25 mocks, there is compatibility within the error associated with the estimation of the standard deviations.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Scatter of the $1\sigma$ statistical errors for the baseline $P+B$ analysis choice on the parameters of interest, $\{\alpha_\parallel, \alpha_\perp, f\sigma_{s8}, m+n\}$ as well as the histogram for the $\chi^2$ distribution normalized by the number of degrees of freedom ($\chi^2/\text{dof}$). The performance for the 25 Abacus AltMTL mocks is displayed in blue circles, whereas the red cross displays the actual DESI DR1 data performance. The orange points with error bars are the errors obtained as the standard deviation of the medians of the 25 individual fits (the error bars corresponding to the error on the standard deviation, estimated as $\sigma/\sqrt{2N-2}$, with $N=25$). The green curves delimit the approximated $1-2\sigma$ regions of the error distribution, estimated from the scatter of the blue data-points. The four main panels represent the redshift bins studied in this paper, as labelled. The reported errors do not contain any systematic error budget, but the EZmock-derived covariance does contain the correcting factors of table 3 of ref. [61]. The errors reported for the DESI data are consistent with those of the Abacus AltMTL mocks, and with the error derived from the medians of the 25 realisations. Only for the LRG3 and QSO redshift bins, the error on the $m+n$ parameter is on the upper side of the distribution of errors, showing a typical $\sim 2\sigma$ offset from the overall distribution. The larger-than-mocks best-fitting $\chi^2$ reported on the data for the LRG2 bin does not have a noticeable effect on the errors of the parameters of interest compared to what is found on the Abacus mocks.</figcaption></figure>
difference in best-fitting $\chi^2$ values does not translate into a difference in the distribution of the errors of the parameters of interest, and therefore it doesn't impact the cosmology results.

We conclude that there is no evidence that the statistical errors are significantly under/over estimated, and that the rescaling factors that we have applied (following table 3 of ref. [61]) are appropriate for our analysis settings. The high values for the minimum $\chi^2$ seen in the data are also seen in the simulations and are possibly due to a combination of effects including imperfections in the model of the bispectrum, imperfections in the calibration of the covariance matrix and possibly the adoption of a Gaussian likelihood as an imperfect approximation. Even if we were to rescale the covariance matrix by a constant factor as to force the minimum $\chi^2$ per dof to be 1, this would not increase the errors by more than ~20%. This is below our minimal threshold to propagate systematics into the final results.

We conclude that the covariance matrix shortcomings we have seen here do not bias the final results nor the reported size of the (statistical + systematic) error-bars in a significant way.

# 6 Conclusions

We have presented the first joint analysis of the power spectrum and bispectrum signals using DESI DR1 data catalogues, focusing on the LRG and QSO samples, and thus probing the evolution of the universe through a redshift range of $0.4 \le z \le 2.1$. The inclusion of the bispectrum allows us to break parameter degeneracies, in particular, to obtain separate constraints for the amplitude of perturbations ($\sigma_{s8}$) and growth rate ($f$) parameters, which are usually constrained only in the combination $f\sigma_{s8}$.

We use a model and framework very similar to that of BOSS and eBOSS analyses [58, 59], which differs from the official DESI pipeline [61] in several aspects:

* We work with the renormalized perturbation theory model (RPT) expansion at two-loop, while [61] uses an EFT-based model at one-loop.

* In the three LRG bins we use the signal from the power spectrum monopole, quadrupole and hexadecapole, while [61] uses only the power spectrum monopole and quadrupole.

* We consider a range of scales of $k < 0.15 h \text{ Mpc}^{-1}$, while [61] considers $k < 0.20 h \text{ Mpc}^{-1}$.

* We do not use all the systematic mitigation strategies that [61] uses, in particular the $\theta$-cut and the marginalization over the non-linearities of systematic weights (although the latter is mostly relevant for ELG galaxies which we do not include).

* We do not use any ELG redshift bin, given that this tracer fails to pass some of our systematic tests. We also do not use the BGS tracer.

* We include the bispectrum monopole.

Apart from these differences, our analysis closely matches that of the official DESI pipeline: we use the same mocks (see section 2.1), both for validation and covariance, the same blinding procedure (section 2), the same approach for the ShapeFit compression (section 3.1), and the same treatment and characterization of systematic errors (section 4).
The remarkable consistency between the P-only results presented here and the official DESI one highlights the exquisite fidelity of DESI data and the astounding robustness of cosmological constraints obtained from compressed-variable analyses of large-scale structure data.

The inclusion of the bispectrum monopole represents the main novelty of this work. We use the GEO-FPT bispectrum model presented in [56], for a range of scales of $0.02 \le k \ [h \ \text{Mpc}^{-1}] \le 0.12$. We limit ourselves to the bispectrum monopole ($B_0$) as we lack of a sufficiently accurate modelling of the effects of the survey window function for the bispectrum multipoles. The adopted bispectrum monopole window convolution is the same as that used in some BOSS and eBOSS analyses [45, 47].

We perform a suite of systematic tests, following closely ref. [61], to determine the size of the systematic error budget and to decide which summary statistics we use in each redshift bin before unblinding the data. The largest sources of systematics come from the modelling and fibre assignment process, with $\sigma_{s8}$ being the most affected parameter in terms of percentage of the DR1 statistical error bar. From the systematic error estimates and the cosmological parameter constraints performed in blinded data, we choose to use as our baseline analysis the combination $P_{024} + B_0$ for the three LRG redshift bins, and only $P_{02} + B_0$ for the QSO tracer.

In this way (and assuming the different redshift bins are uncorrelated), we obtain cumulative constraints of 10.1% and 8.4% for respectively $f$ and $\sigma_{s8}$. Our constraints for the parameters $\alpha_{\text{iso}}$, $\alpha_{\text{AP}}$, $f\sigma_{s8}$, $m+n$ are consistent (and competitive) with the ones obtained in the main DESI collaboration ShapeFit analysis [61]. In particular, the inclusion of the bispectrum monopole reduces the error bars of $\alpha_{\text{iso}}$, $f\sigma_{s8}$, $m+n$ by respectively 9%, 9% and 11% with respect to our power spectrum-only analysis.

We envision this to be a first milestone towards the promise of tight cosmological constraints (including $f$ and $\sigma_{s8}$ when considered separately) with the combination of two-point and higher-order statistics. In particular, the ability to consider parameter constraints on $f$ and $\sigma_{s8}$ separately will allow us to constrain modifications of general relativity (where the relationship $f = \Omega_m^\gamma$ with $\gamma \approx 6/11$ no longer holds [117, 118]). Aside from the obvious improvement in precision that will be possible with the upcoming DESI data releases, improvements will include strategies to mitigate the fibre assignment systematic errors and the inclusion of the bispectrum quadrupoles in the analysis. This is a major challenge since the window function treatment is highly non-trivial, but the extra anisotropic information present may provide significant enhancement of the precision and accuracy in most cosmological parameters, as we saw in [56].

In a forthcoming work [79], we will provide the interpretation of the compressed parameter constraints presented here in terms of traditional cosmological parameters for the $\Lambda$CDM model and its extensions.

## Acknowledgments

SNM, HGM and LV thank Pauline Zarrouk, Caroline Guandalin, Alex Krolewski and Ruiyang Zhao for the helpful discussion and suggestions. SNM acknowledges funding from the official doctoral programme of the University of Barcelona for the development of a research project under the PREDOCS-UB grant. HGM acknowledges support through the Leonardo programme (LEO23-1-897) of the BBVA foundation and through the programmes Ramón y
Cajal (RYC-2021-034104) and Consolidación Investigadora (CNS2023-144605) of the Spanish Ministry of Science and Innovation.

Funding for this work was partially provided by the Spanish MINECO under project PID2022-141125NB-I00 MCIN/AEI, and the “Center of Excellence Maria de Maeztu 2020-2023” award to the ICCUB (CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033).

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies.

The authors are honored to be permitted to conduct scientific research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.

This work has made use of the following publicly available codes: GEO-FPT [56], RUSTICO [59], BRASS [59], EMCEE [119], GSL [120], SciPy [121], NumPy [122], GETDIST [123], ASTROPY [124], MATPLOTLIB [125]. We are grateful to the developers who made these codes public.

## A Perturbation theory modelling of the power spectrum and bispectrum

In this section, we review the perturbation theory theoretical formalism underlying the modelling used here for the power spectrum and bispectrum.

The redshift space galaxy power spectrum, $P_g$, is computed from the non-linear matter power spectrum, $P_{g,\delta\delta}$ (also referred to as $P_{\text{NL}}$ in the main text and in refs. [56, 126]), the density-velocity, $P_{g,\delta\theta}$, and velocity-velocity, $P_{g,\theta\theta}$, power spectra, according to the TNS model [127, 128],

$$P_g(k, \mu) = D_{\text{FoG}}^P(k, \mu, \sigma_P) [P_{g,\delta\delta}(k) + 2f\mu^2 P_{g,\delta\theta}(k) + f^2\mu^4 P_{\theta\theta}(k) \\ + b_1^2 A^{\text{TNS}}(k, \mu, f/b_1) + b_1^4 B^{\text{TNS}}(k, \mu, f/b_1)], \quad (A.1)$$

where $f$ denotes the logarithmic growth rate of perturbations, and $d\ln\delta/d\ln a$ and $P_{g,\delta\delta}, P_{g,\delta\theta}$ are computed as in [129] using the RPT terms of [57] for $P_{\delta\delta}, P_{\delta\theta}, P_{\theta\theta}$. In doing so, we are using the bias expansion {$b_1, b_2, b_{s^2}, b_{3\text{nl}}$} and assuming the Lagrangian local bias relations [92, 93, 129, 130],

$$b_{s^2} = -\frac{4}{7}(b_1 - 1); \quad b_{3\text{nl}} = \frac{32}{315}(b_1 - 1). \qquad (A.2)$$
Additionally, the functions $A^{\text{TNS}}$, $B^{\text{TNS}}$ are defined in [127], $\mu$ is the cosine of the angle of $k$ with the line of sight, and $D_{\text{FoG}}^{P}$ is a damping factor that accounts for the Fingers-of-God (FoG) effect of redshift space distortions (RSD) [131]. We model the FoG damping factor for the power spectrum as,

$$D_{\text{FoG}}^{P}(k, \mu, \sigma_{\text{FoG}}^{P}) = \frac{1}{(1 + k^2 \mu^2 \sigma_{\text{P}}^2 / 2)^2}, \qquad (\text{A.3})$$

where $\sigma_P$ is a free parameter to be constrained by the data.

The tree-level redshift space bispectrum can then be written in the following way at tree-level order:

$$B^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = D_{\text{FoG}}^{B}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) [2Z_{1}^{\text{SPT}}(\mathbf{k}_1)Z_{2}^{\text{SPT}}(\mathbf{k}_2)Z_{3}^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2)P_{\text{L}}(k_1)P_{\text{L}}(k_2) + 2\text{perm.}], \quad (\text{A.4})$$

where the kernels $Z_1^{\text{SPT}}$, $Z_2^{\text{SPT}}$ are computed as [25, 132, 133],

$$Z_1^{\text{SPT}}(\mathbf{k}) = b_1 + f\mu^2, \\ Z_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) = b_1 F_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) + f\mu_{12}^2 G_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) + \frac{b_1 f}{2} \left( \mu_1^2 + \mu_2^2 + \mu_1 \mu_2 \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) \right) \\ + f^2 \mu_1 \mu_2 \left( \mu_1 \mu_2 + \frac{1}{2} \left( \mu_1^2 \frac{k_1}{k_2} + \mu_2^2 \frac{k_2}{k_1} \right) \right) + \frac{1}{2} \left( b_2 + b_{s2} S_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) \right), \quad (\text{A.5})$$

where $\mu_{ij} \equiv (k_i \mu_i + k_j \mu_j)/|\mathbf{k}_i + \mathbf{k}_j|$. The $G_2^{\text{SPT}}$ and $S_2^{\text{SPT}}$ kernels in Standard Perturbation Theory (SPT) are given by,

$$G_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) = \frac{3}{7} + \frac{1}{2} \cos(\theta_{12}) \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) + \frac{4}{7} \cos^2(\theta_{12}), \quad (\text{A.6})$$

$$S_2^{\text{SPT}}(\mathbf{k}_1, \mathbf{k}_2) = \cos(\theta_{12})^2 - \frac{1}{3}. \qquad (\text{A.7})$$

Additionally, for the bispectrum, we parametrise the FoG damping factor as [25, 102],

$$D_{\text{FoG}}^{B}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = (1 + [k_1^2 \mu_1^2 + k_2^2 \mu_2^2 + k_3^2 \mu_3^2]^2 \sigma_B^4 / 2)^{-2}, \quad (\text{A.8})$$

again with $\sigma_B$ being a free parameter.

We model the deviations from Poissonian shot-noise with the parameters $A_P$, $A_B$, which modify the Poisson prediction as in [47, 134],

$$P_{\text{noise}} = \left( 1 - \frac{A_P}{\alpha_{\parallel} \alpha_{\perp}^2} \right) P_{\text{Poisson}}, \quad (\text{A.9})$$

$$B_{\text{noise}}(k_1, k_2, k_3) = \left( 1 - \frac{A_B}{\alpha_{\parallel}^2 \alpha_{\perp}^4} \right) B_{\text{Poisson}}(k_1, k_2, k_3). \quad (\text{A.10})$$

The power spectrum and bispectrum redshift space multipoles are then obtained by integrating the expansion of the power spectrum and bispectrum dependence on the angle
with respect to the line of sight in terms of Legendre polynomials $L_i$, so that

$$P_{\ell}(k) = \frac{2\ell + 1}{2\alpha_{\parallel}\alpha_{\perp}^2} \int_{-1}^{1} d\mu P(k, \mu)L_{\ell}(\mu), \qquad (A.11)$$

$$B_{\ell_i}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \frac{2\ell + 1}{4\pi\alpha_{\parallel}^2\alpha_{\perp}^4} \int_{-1}^{1} d\mu_1 \int_{0}^{2\pi} d\phi B(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3)L_{\ell}(\mu_i). \qquad (A.12)$$

Here $\phi$ is defined as the angle fulfilling $\mu_2 = \mu_1 \cos \theta_{12} - \sqrt{(1 - \mu_1^2)(1 - \cos \theta_{12}^2)} \cos \phi$, and $\ell_i$ refers to the multipole of order $\ell$ ($\ell = 0, 2$ corresponding respectively to the monopole and quadrupole). For the bispectrum quadrupoles ($\ell = 2$), $\ell_i$ is the quadrupole corresponding to integrating over the Legendre polynomial applied to the cosine of the $i$-th angle. In this work, we only use the bispectrum monopole, due to the current limitations of the modelling of the window function. The power spectrum multipole expansion of (A.11) was proposed in [135, 136], while the bispectrum expansion and choice of variables was first used in [102].

## B Blinded data

The blinded cosmological constraints for the parameters $\{\alpha_{\text{iso}}, \alpha_{\text{AP}}, f\sigma_{\text{s8}}, m + n, f, \sigma_{\text{s8}}\}$ for the four combinations of summary statistics that we consider in this work ($P_{02}, P_{024}, P_{02} + B_0, P_{024} + B_0$) are shown together in figure 13.

These results for the blinded data determine the combination of summary statistics to use after unblinding. We aim at striking a balance between two aspects:

* **Statistical errors:** In general, adding more multipoles or the bispectrum to the data-vector decreases the error bars

* **Systematic errors:** Some scales or summary statistics are more prone to systematic errors than others. The systematic errors presented in tables 4, 5, 9, 10 allow us to quantify this.

We wish to select data-vectors which add significant information (and thus reduce the statistical error bars) while not increasing the systematic error budget.

Figure 13 displays the constraints on the ShapeFit cosmological parameters for the different DESI redshift bins, and the different combinations of summary statistics considered (symbols in different colours). The addition of the power spectrum hexadecapole, $P_4$, significantly enhances the LRG constraints for the anisotropic parameter $\alpha_{\text{AP}}$, and, to a lesser extent, $f\sigma_{\text{s8}}$. The comparison of tables 9 vs 10 and tables 4 vs 5 shows that for LRGs the inclusion of $P_4$ does not increase the systematic error budget in any significant way.

Conversely, in the case of the QSO tracers, the inclusion of $P_4$ does not reduce appreciably the statistical errors but increases notably the systematic error budget (specifically for the parameters $f\sigma_{\text{s8}}$ and $\alpha_{\text{AP}}$). This is not unexpected, since the QSO tracers may not have sufficient number density for the signal-to-noise of the power spectrum hexadecapole to be big enough.

The inclusion of the bispectrum, aside from breaking the $f\sigma_{\text{s8}}$ degeneracy, reduces the error bars of the $\alpha_{\text{iso}}$ and $m + n$ parameters. Additionally, the $P_{024} + B_0$ multipoles
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Cosmological constraints obtained from the four blinded redshift bins (LRG1, LRG2, LRG3, QSO), by using each of the four combinations of summary statistics that we consider in this work: $P_{02}$, $P_{024}$, $P_{02} + B_0$, $P_{024} + B_0$. The parameters $f$ and $\sigma_{s8}$ are not shown separately, neither for the power spectrum, as usual, nor for the QSO redshift bin, where we saw that even in an analysis including the bispectrum the two parameters were markedly degenerate.</figcaption></figure>
<table><caption>**Table 7.** Results for the nuisance parameters from our two main analyses, $P+B$ and $P$, for the LRG and QSO redshift bins, analogue to table 6. Each row displays the maximum a posteriori (MAP), together with the $1\sigma$ region (which in this case only accounts for the statistical errors, as there is no systematic contribution calculated) centred in the MAP value. In the cases with significant projection effects, the MAP may fall outside of the $1\sigma$ region. We additionally report the corresponding values for $\chi^2$ over the degrees of freedom for each case.</caption><thead><tr><th></th><th>Sample</th><th>b<sub>1</sub>σ<sub>s8</sub></th><th>b<sub>2</sub>σ<sub>s8</sub><sup>3</sup></th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>A<sub>P</sub></th><th>A<sub>B</sub></th><th>σ<sub>P</sub> [Mpc h<sup>-1</sup>]</th><th>σ<sub>B</sub> [Mpc h<sup>-1</sup>]</th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td rowspan="4">P + B</td><td>LRG1</td><td>1.159<sup>+0.084</sup><sub>-0.024</sub></td><td>0.213<sup>+0.090</sup><sub>-0.044</sub></td><td>1.525<sup>+0.180</sup><sub>-0.074</sub></td><td>0.484<sup>+0.356</sup><sub>-0.129</sub></td><td>0.970<sup>+0.161</sup><sub>-0.082</sub></td><td>0.884<sup>+0.265</sup><sub>-0.102</sub></td><td>4.404<sup>+0.497</sup><sub>-0.386</sub></td><td>0.521<sup>+2.506</sup><sub>-1.010</sub></td><td>241/178</td></tr><tr><td>LRG2</td><td>1.140<sup>+0.040</sup><sub>-0.031</sub></td><td>0.183<sup>+0.058</sup><sub>-0.069</sub></td><td>1.882<sup>+0.169</sup><sub>-0.274</sub></td><td>0.822<sup>+0.487</sup><sub>-0.489</sub></td><td>0.724<sup>+0.203</sup><sub>+0.003</sub></td><td>0.787<sup>+0.317</sup><sub>-0.066</sub></td><td>4.253<sup>+0.335</sup><sub>-0.427</sub></td><td>3.672<sup>+1.364</sup><sub>-2.598</sub></td><td>291/178</td></tr><tr><td>LRG3</td><td>1.133<sup>+0.030</sup><sub>-0.026</sub></td><td>0.233<sup>+0.060</sup><sub>-0.014</sub></td><td>2.265<sup>+0.551</sup><sub>-0.106</sub></td><td>1.863<sup>+2.251</sup><sub>-0.197</sub></td><td>0.863<sup>+0.155</sup><sub>+0.001</sub></td><td>0.998<sup>+0.153</sup><sub>-0.102</sub></td><td>3.339<sup>+0.444</sup><sub>-0.355</sub></td><td>4.511<sup>+0.163</sup><sub>-3.609</sub></td><td>253/178</td></tr><tr><td>QSO</td><td>0.890<sup>+0.069</sup><sub>-0.007</sub></td><td>0.058<sup>+0.052</sup><sub>-0.011</sub></td><td>—</td><td>—</td><td>0.981<sup>+0.118</sup><sub>-0.013</sub></td><td>1.015<sup>+0.231</sup><sub>-0.040</sub></td><td>4.563<sup>+0.769</sup><sub>-0.340</sub></td><td>7.526<sup>+1.990</sup><sub>-5.530</sub></td><td>193/166</td></tr><tr><td rowspan="4">P</td><td>LRG1</td><td>1.136<sup>+0.002</sup><sub>-0.078</sub></td><td>0.260<sup>+0.500</sup><sub>-0.400</sub></td><td>—</td><td>—</td><td>0.772<sup>+0.427</sup><sub>+0.118</sub></td><td>—</td><td>4.057<sup>-0.263</sup><sub>-1.192</sub></td><td>—</td><td>33/31</td></tr><tr><td>LRG2</td><td>1.112<sup>+0.023</sup><sub>-0.054</sub></td><td>0.621<sup>+0.139</sup><sub>-0.761</sub></td><td>—</td><td>—</td><td>0.696<sup>+0.503</sup><sub>+0.194</sub></td><td>—</td><td>4.083<sup>-0.289</sup><sub>-1.218</sub></td><td>—</td><td>37/31</td></tr><tr><td>LRG3</td><td>1.113<sup>+0.021</sup><sub>-0.055</sub></td><td>0.429<sup>+0.331</sup><sub>-0.569</sub></td><td>—</td><td>—</td><td>0.911<sup>+0.288</sup><sub>-0.021</sub></td><td>—</td><td>3.258<sup>+0.536</sup><sub>-0.392</sub></td><td>—</td><td>31/31</td></tr><tr><td>QSO</td><td>0.904<sup>+0.026</sup><sub>-0.043</sub></td><td>0.047<sup>+0.293</sup><sub>-0.213</sub></td><td>—</td><td>—</td><td>0.992<sup>+0.115</sup><sub>-0.031</sub></td><td>—</td><td>4.731<sup>+0.973</sup><sub>-0.355</sub></td><td>—</td><td>20/18</td></tr></tbody></table>

combination features slightly less projection effects than the case without the power spectrum hexadecapole, particularly in the $\alpha_{AP}$, $m+n$, $f$ and $\sigma_{s8}$ parameters.<fnref n="21" />

Taking all these points into account, we define our baseline set of multipoles as $P_{024} + B_0$ for the three LRG redshift bins, and $P_{02} + B_0$ for the QSOs. Our baseline power spectrum-only analysis will likewise consist in $P_{024}$ for the LRGs and $P_{02}$ for QSOs.

## C Nuisance parameters

In this section, we report the MAP values for the nuisance parameters obtained from the analysis described in section 5. Table 7 displays the results for the galaxy bias parameters, the shot noise amplitudes and the Fingers-of-God parameters for both types of analysis, with and without the bispectrum signal, as previously done in table 6 for the ShapeFit cosmological parameters. The baseline choices of those results are the same to those presented before in section 5: for the LRG samples we vary $\sigma_{s8}$ and choose to report both $\{b_1, b_2\}$, as well as the combination of parameters best-measured, $\{b_1\sigma_{s8}, b_2\sigma_{s8}^3\}$. On the other hand, for the QSO and for the LRG samples when only $P$ is employed, we keep $\sigma_{s8}$ fixed, and re-interpret the best-fitting bias values as a product of $b_1$ and a power of $\sigma_{s8}$. We determine that the combination that for the power spectrum and bispectrum that keeps $b_x\sigma_{s8}^n$ uncorrelated with $\sigma_{s8}$ is $n=1$ for $b_1$ (as expected from the Kaiser limit at large scales), and $n=3$ for $b_2$, as a resulting of combining the power spectrum and bispectrum. We then choose those combination of variables to be reported here.

Figure 14 displays these same results in a triangle-plot style.

<footnote n="21">²¹The presence of projection effects can be inferred from the displacement of the maximum a posteriori (MAP) value from the centre of the error bars. When this shift is significant, it indicates that the data provide weak constraints on the parameter, since the posterior distribution is skewed.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Posteriors obtained for the ShapeFit analysis of $P+B$ for the nuisance parameters $b_1$, $b_2$, $A_P$, $A_B$, $\sigma_P$ and $\sigma_B$ (units of $\sigma_{P,B}$ are Mpch$^{-1}$), for each sample as labelled. The cosmological parameters are also varied in this fit (see figure 11), but not shown for clarity. We also report the galaxy bias results in terms of $b_1\sigma_{s8}$ and $b_2\sigma_{s8}^3$, as these are the effective combinations better constrained for the $P+B$ case. The MAPs values of these contours are reported in table 7.</figcaption></figure>

D Covariances and additional tables

We will report here, upon acceptance, the covariance of the cosmological parameters for each tracer.

We also report here some numerical details not included in the main text. The cosmological parameter constraints for the non-baseline choices of data-vector (i.e. the ones involving $P_{024}$ in QSOs and $P_{02}$ in LRGs, for both power spectrum and joint power spectrum-bispectrum analyses), are shown in table 8.
<table><caption>**Table 8.** Analogous to table 6, in this case showing the combinations of statistics that are not a part of our baseline analysis.</caption><thead><tr><th></th><th>Sample</th><th>α<sub>iso</sub></th><th>α<sub>AP</sub></th><th>dfσ<sub>s8</sub></th><th>m + n</th><th>df</th><th>dσ<sub>s8</sub></th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td rowspan="3">P<sub>02</sub> + B</td><td>LRG1</td><td>0.998<sup>+0.018</sup><sub>-0.020</sub></td><td>0.961<sup>+0.106</sup><sub>-0.036</sub></td><td>1.137<sup>+0.131</sup><sub>-0.119</sub></td><td>0.001<sup>+0.020</sup><sub>-0.0795</sub></td><td>1.014<sup>+0.129</sup><sub>-0.217</sub></td><td>1.121<sup>+0.197</sup><sub>-0.049</sub></td><td>231/165</td></tr><tr><td>LRG2</td><td>0.967<sup>+0.029</sup><sub>-0.014</sub></td><td>1.064<sup>+0.136</sup><sub>-0.093</sub></td><td>0.993<sup>+0.134</sup><sub>-0.164</sub></td><td>0.035<sup>+0.023</sup><sub>-0.072</sub></td><td>0.745<sup>+0.379</sup><sub>-0.173</sub></td><td>1.332<sup>+0.200</sup><sub>-0.431</sub></td><td>271/165</td></tr><tr><td>LRG3</td><td>1.002<sup>+0.022</sup><sub>-0.014</sub></td><td>1.093<sup>+0.085</sup><sub>-0.083</sub></td><td>0.896<sup>+0.153</sup><sub>-0.111</sub></td><td>0.036<sup>+0.009</sup><sub>-0.085</sub></td><td>0.699<sup>+0.418</sup><sub>-0.024</sub></td><td>1.281<sup>-0.042</sup><sub>-0.404</sub></td><td>240/165</td></tr><tr><td>P<sub>024</sub> + B</td><td>QSO</td><td>1.019<sup>+0.026</sup><sub>-0.022</sub></td><td>1.082<sup>+0.029</sup><sub>-0.053</sub></td><td>1.016<sup>+0.093</sup><sub>-0.091</sub></td><td>0.067<sup>+0.036</sup><sub>-0.046</sub></td><td>—</td><td>—</td><td>213/179</td></tr><tr><td rowspan="3">P<sub>02</sub></td><td>LRG1</td><td>0.965<sup>+0.026</sup><sub>-0.016</sub></td><td>1.042<sup>+0.085</sup><sub>-0.075</sub></td><td>0.935<sup>+0.206</sup><sub>-0.047</sub></td><td>0.009<sup>+0.010</sup><sub>-0.097</sub></td><td>—</td><td>—</td><td>24/18</td></tr><tr><td>LRG2</td><td>0.949<sup>+0.022</sup><sub>-0.015</sub></td><td>1.035<sup>+0.078</sup><sub>-0.072</sub></td><td>0.993<sup>+0.164</sup><sub>-0.089</sub></td><td>0.016<sup>+0.032</sup><sub>-0.076</sub></td><td>—</td><td>—</td><td>20/18</td></tr><tr><td>LRG3</td><td>0.994<sup>+0.020</sup><sub>-0.019</sub></td><td>1.056<sup>+0.096</sup><sub>-0.070</sub></td><td>0.951<sup>+0.172</sup><sub>-0.127</sub></td><td>-0.011<sup>+0.037</sup><sub>-0.067</sub></td><td>—</td><td>—</td><td>17/18</td></tr><tr><td>P<sub>024</sub></td><td>QSO</td><td>1.006<sup>+0.031</sup><sub>-0.022</sub></td><td>1.082<sup>+0.029</sup><sub>-0.053</sub></td><td>1.016<sup>+0.093</sup><sub>-0.091</sub></td><td>0.067<sup>+0.036</sup><sub>-0.046</sub></td><td>—</td><td>—</td><td>33/31</td></tr></tbody></table>

<table><caption>**Table 9.** Same as table 4, with the only difference that the summary statistics considered here are the power spectrum monopole, quadrupole and hexadecapole (P<sub>024</sub>), without the bispectrum.</caption><thead><tr><th>[%σ<sub>DR1</sub>]</th><th>P<sub>024</sub></th><th>σ<sub>α<sub>iso</sub></sub></th><th>σ<sub>α<sub>AP</sub></sub></th><th>σ<sub>fσ<sub>s8</sub></sub></th><th>σ<sub>m+n</sub></th></tr></thead><tbody><tr><td rowspan="4">Modelling</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 20</td></tr><tr><td>LRG2</td><td>~ 30</td><td>~ 20</td><td>~ 30</td><td>~ 20</td></tr><tr><td>LRG3</td><td>~ 40</td><td>&lt; 20</td><td>~ 20</td><td>~ 40</td></tr><tr><td>QSO</td><td>~ 40</td><td>&lt; 10</td><td>~ 20</td><td>~ 20</td></tr><tr><td rowspan="4">Observational</td><td>LRG1</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 10</td></tr><tr><td>LRG2</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 10</td></tr><tr><td>LRG3</td><td>&lt; 10</td><td>~ 60</td><td>~ 80</td><td>&lt; 10</td></tr><tr><td>QSO</td><td>~ 50</td><td>~ 50</td><td>~ 130</td><td>&lt; 20</td></tr><tr><td rowspan="4">HOD</td><td>LRG1</td><td>&lt; 10(11)</td><td>&lt; 10(7)</td><td>&lt; 10(9)</td><td>&lt; 10(9)</td></tr><tr><td>LRG2</td><td>&lt; 10(16)</td><td>&lt; 10(9)</td><td>&lt; 10(14)</td><td>&lt; 20(13)</td></tr><tr><td>LRG3</td><td>&lt; 10(16)</td><td>&lt; 10(12)</td><td>&lt; 10(17)</td><td>&lt; 10(14)</td></tr><tr><td>QSO</td><td>&lt; 10(4)</td><td>&lt; 10(7)</td><td>&lt; 10(15)</td><td>&lt; 20(29)</td></tr><tr><td rowspan="4">Fiducial cosmology</td><td>LRG1</td><td>&lt; 10(17)</td><td>&lt; 10(7)</td><td>&lt; 10(13)</td><td>&lt; 20(45)</td></tr><tr><td>LRG2</td><td>&lt; 10(25)</td><td>&lt; 10(9)</td><td>&lt; 20(19)</td><td>&lt; 20(61)</td></tr><tr><td>LRG3</td><td>&lt; 10(25)</td><td>&lt; 10(12)</td><td>&lt; 20(23)</td><td>&lt; 20(65)</td></tr><tr><td>QSO</td><td>&lt; 10(24)</td><td>&lt; 10(10)</td><td>&lt; 20(40)</td><td>&lt; 20(55)</td></tr><tr><td rowspan="4">Total</td><td>LRG1</td><td>&lt; 20</td><td>&lt; 20</td><td>~ 20</td><td>&lt; 20</td></tr><tr><td>LRG2</td><td>~ 30</td><td>~ 20</td><td>~ 50</td><td>~ 20</td></tr><tr><td>LRG3</td><td>~ 40</td><td>~ 60</td><td>~ 90</td><td>~ 40</td></tr><tr><td>QSO</td><td>~ 64</td><td>~ 50</td><td>~ 132</td><td>~ 20</td></tr></tbody></table>

The systematic error tables for the power spectrum analysis (not including the bispectrum), analogous to tables 4 and 5, are respectively table 9 and 10.
<table><caption>**Table 10.** Same as table 4, with the only difference that the summary statistics considered here are the power spectrum monopole and quadrupole ($P_{02}$) only.</caption><thead><tr><th>[%σ<sub>DR1</sub>]</th><th>P<sub>02</sub></th><th>σ<sub>α<sub>iso</sub></sub></th><th>σ<sub>α<sub>AP</sub></sub></th><th>σ<sub>fσ<sub>s8</sub></sub></th><th>σ<sub>m+n</sub></th></tr></thead><tbody><tr><td rowspan="4">Modelling</td><td>LRG1</td><td>~ 20</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 10</td></tr><tr><td>LRG2</td><td>~ 30</td><td>&lt; 10</td><td>&lt; 20</td><td>&lt; 20</td></tr><tr><td>LRG3</td><td>~ 40</td><td>&lt; 20</td><td>~ 40</td><td>~ 40</td></tr><tr><td>QSO</td><td>~ 50</td><td>&lt; 20</td><td>&lt; 10</td><td>&lt; 20</td></tr><tr><td rowspan="4">Observational</td><td>LRG1</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 10</td></tr><tr><td>LRG2</td><td>&lt; 10</td><td>&lt; 20</td><td>&lt; 20</td><td>&lt; 10</td></tr><tr><td>LRG3</td><td>&lt; 10</td><td>&lt; 10</td><td>&lt; 20</td><td>&lt; 10</td></tr><tr><td>QSO</td><td>~ 50</td><td>&lt; 10</td><td>~ 70</td><td>&lt; 20</td></tr><tr><td rowspan="4">HOD</td><td>LRG1</td><td>&lt; 10(10)</td><td>&lt; 10(6)</td><td>&lt; 10(8)</td><td>&lt; 10(13)</td></tr><tr><td>LRG2</td><td>&lt; 10(15)</td><td>&lt; 10(9)</td><td>&lt; 10(12)</td><td>&lt; 10(18)</td></tr><tr><td>LRG3</td><td>&lt; 10(15)</td><td>&lt; 10(10)</td><td>&lt; 10(13)</td><td>&lt; 10(20)</td></tr><tr><td>QSO</td><td>&lt; 10(4)</td><td>&lt; 10(8)</td><td>&lt; 10(10)</td><td>&lt; 20(27)</td></tr><tr><td rowspan="4">Fiducial cosmology</td><td>LRG1</td><td>&lt; 10(18)</td><td>&lt; 10(7)</td><td>&lt; 10(15)</td><td>&lt; 20(43)</td></tr><tr><td>LRG2</td><td>&lt; 10(28)</td><td>&lt; 10(10)</td><td>&lt; 20(22)</td><td>~ 20(59)</td></tr><tr><td>LRG3</td><td>&lt; 10(26)</td><td>&lt; 10(11)</td><td>&lt; 20(25)</td><td>~ 20(65)</td></tr><tr><td>QSO</td><td>&lt; 20(28)</td><td>&lt; 10(6)</td><td>&lt; 20(36)</td><td>&lt; 20(54)</td></tr><tr><td rowspan="4">Total</td><td>LRG1</td><td>~ 20</td><td>&lt; 20</td><td>~ 30</td><td>&lt; 20</td></tr><tr><td>LRG2</td><td>~ 30</td><td>&lt; 20</td><td>~ 40</td><td>~ 20</td></tr><tr><td>LRG3</td><td>~ 40</td><td>&lt; 20</td><td>~ 57</td><td>~ 45</td></tr><tr><td>QSO</td><td>~ 71</td><td>&lt; 20</td><td>~ 70</td><td>&lt; 20</td></tr></tbody></table>

**Data Availability Statement.** This article has associated data in a data repository. All data from the tables and figures will be available in machine-readable format at https://zenodo.org/records/14944381 upon acceptance in compliance with the DESI data management plan.

References

[1] L. Verde, T. Treu and A.G. Riess, *Tensions between the early and late universe*, Nature Astron. 3 (2019) 891 [arXiv:1907.10625] [INSPIRE].

[2] SDSS collaboration, *The 2.5 m Telescope of the Sloan Digital Sky Survey*, Astron. J. 131 (2006) 2332 [astro-ph/0602326] [INSPIRE].

[3] 2DFGRS collaboration, *The 2dF Galaxy Redshift Survey: The power spectrum and the matter content of the Universe*, Mon. Not. Roy. Astron. Soc. 327 (2001) 1297 [astro-ph/0105252] [INSPIRE].

[4] BOSS collaboration, *The Baryon Oscillation Spectroscopic Survey of SDSS-III*, Astron. J. 145 (2013) 10 [arXiv:1208.0022] [INSPIRE].

[5] EBOSS collaboration, *The SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Overview and Early Data*, Astron. J. 151 (2016) 44 [arXiv:1508.04473] [INSPIRE].

[6] DESI collaboration, *The DESI Experiment*, a whitepaper for Snowmass 2013, arXiv:1308.0847 [INSPIRE].
[7] DESI collaboration, *The DESI Experiment Part I: Science, Targeting, and Survey Design*, arXiv:1611.00036 [INSPIRE].

[8] DESI collaboration, *The DESI Experiment Part II: Instrument Design*, arXiv:1611.00037 [INSPIRE].

[9] DESI collaboration, *Overview of the Instrumentation for the Dark Energy Spectroscopic Instrument*, Astron. J. 164 (2022) 207 [arXiv:2205.10939] [INSPIRE].

[10] DESI collaboration, *The Robotic Multiobject Focal Plane System of the Dark Energy Spectroscopic Instrument (DESI)*, Astron. J. 165 (2023) 9 [arXiv:2205.09014] [INSPIRE].

[11] DESI collaboration, *The Optical Corrector for the Dark Energy Spectroscopic Instrument*, Astron. J. 168 (2024) 95 [arXiv:2306.06310] [INSPIRE].

[12] DESI collaboration, *The Spectroscopic Data Processing Pipeline for the Dark Energy Spectroscopic Instrument*, Astron. J. 165 (2023) 144 [arXiv:2209.14482] [INSPIRE].

[13] DESI collaboration, *Survey Operations for the Dark Energy Spectroscopic Instrument*, Astron. J. 166 (2023) 259 [arXiv:2306.06309] [INSPIRE].

[14] DESI collaboration, *Target Selection and Validation of DESI Luminous Red Galaxies*, Astron. J. 165 (2023) 58 [arXiv:2208.08515] [INSPIRE].

[15] DESI collaboration, *Validation of the Scientific Program for the Dark Energy Spectroscopic Instrument*, Astron. J. 167 (2024) 62 [arXiv:2306.06307] [INSPIRE].

[16] DESI collaboration, *The Early Data Release of the Dark Energy Spectroscopic Instrument*, Astron. J. 168 (2024) 58 [arXiv:2306.06308] [INSPIRE].

[17] DESI collaboration, *DESI 2024 III: baryon acoustic oscillations from galaxies and quasars*, JCAP 04 (2025) 012 [arXiv:2404.03000] [INSPIRE].

[18] DESI collaboration, *DESI 2024 IV: Baryon Acoustic Oscillations from the Lyman alpha forest*, JCAP 01 (2025) 124 [arXiv:2404.03001] [INSPIRE].

[19] DESI collaboration, *DESI 2024 II: Sample Definitions, Characteristics, and Two-point Clustering Statistics*, arXiv:2411.12020 [INSPIRE].

[20] DESI collaboration, *DESI 2024 VI: cosmological constraints from the measurements of baryon acoustic oscillations*, JCAP 02 (2025) 021 [arXiv:2404.03002] [INSPIRE].

[21] C. Poppett et al., *Overview of the Fiber System for the Dark Energy Spectroscopic Instrument*, Astron. J. 168 (2024) 245.

[22] T. Baldauf, U. Seljak and L. Senatore, *Primordial non-Gaussianity in the Bispectrum of the Halo Density Field*, JCAP 04 (2011) 006 [arXiv:1011.1513] [INSPIRE].

[23] A. Moradinezhad Dizgah et al., *Primordial Non-Gaussianity from Biased Tracers: Likelihood Analysis of Real-Space Power Spectrum and Bispectrum*, JCAP 05 (2021) 015 [arXiv:2010.14523] [INSPIRE].

[24] D. Gualdi, S. Novell-Masot, H. Gil-Marín and L. Verde, *Matter trispectrum: theoretical modelling and comparison to N-body simulations*, JCAP 01 (2021) 015 [arXiv:2009.02290] [INSPIRE].

[25] L. Verde, A.F. Heavens, S. Matarrese and L. Moscardini, *Large scale bias in the universe. 2. Redshift space bispectrum*, Mon. Not. Roy. Astron. Soc. 300 (1998) 747 [astro-ph/9806028] [INSPIRE].
[26] S. Matarrese, L. Verde and A.F. Heavens, *Large-scale bias in the Universe: Bispectrum method*, Mon. Not. Roy. Astron. Soc. **290** (1997) 651 [astro-ph/9706059] [INSPIRE].

[27] E. Sefusatti, M. Crocce, S. Pueblas and R. Scoccimarro, *Cosmology and the Bispectrum*, Phys. Rev. D **74** (2006) 023522 [astro-ph/0604505] [INSPIRE].

[28] J.N. Fry and E. Gaztañaga, *Biasing and hierarchical statistics in large-scale structure*, Astrophys. J. **413** (1993) 447 [astro-ph/9302009] [INSPIRE].

[29] R. Ruggeri, E. Castorina, C. Carbone and E. Sefusatti, *DEMUni: Massive neutrinos and the bispectrum of large scale structures*, JCAP **03** (2018) 003 [arXiv:1712.02334] [INSPIRE].

[30] C.H. Hahn, F. Villaescusa-Navarro, E. Castorina and R. Scoccimarro, *Constraining $M_\nu$ with the bispectrum. Part I. Breaking parameter degeneracies*, JCAP **03** (2020) 040 [arXiv:1909.11107] [INSPIRE].

[31] A. Oddo et al., *Cosmological parameters from the likelihood analysis of the galaxy power spectrum and bispectrum in real space*, JCAP **11** (2021) 038 [arXiv:2108.03204] [INSPIRE].

[32] N. Bartolo, E. Bellini, D. Bertacca and S. Matarrese, *Matter bispectrum in cubic Galileon cosmologies*, JCAP **03** (2013) 034 [arXiv:1301.4831] [INSPIRE].

[33] E. Bellini, R. Jimenez and L. Verde, *Signatures of Horndeski gravity on the Dark Matter Bispectrum*, JCAP **05** (2015) 057 [arXiv:1504.04341] [INSPIRE].

[34] D. Bertacca et al., *Relativistic wide-angle galaxy bispectrum on the light-cone*, Phys. Rev. D **97** (2018) 023531 [arXiv:1705.09306] [INSPIRE].

[35] H. Gil-Marín et al., *The Bispectrum of $f(R)$ Cosmologies*, JCAP **11** (2011) 019 [arXiv:1109.2115] [INSPIRE].

[36] V. Yankelevich and C. Porciani, *Cosmological information in the redshift-space bispectrum*, Mon. Not. Roy. Astron. Soc. **483** (2019) 2078 [arXiv:1807.07076] [INSPIRE].

[37] W.R. Coulton et al., *Constraining Neutrino Mass with the Tomographic Weak Lensing Bispectrum*, JCAP **05** (2019) 043 [arXiv:1810.02374] [INSPIRE].

[38] P. Gagrani and L. Samushia, *Information Content of the Angular Multipoles of Redshift-Space Galaxy Bispectrum*, Mon. Not. Roy. Astron. Soc. **467** (2017) 928 [arXiv:1610.03488] [INSPIRE].

[39] D. Gualdi and L. Verde, *Galaxy redshift-space bispectrum: the Importance of Being Anisotropic*, JCAP **06** (2020) 041 [arXiv:2003.12075] [INSPIRE].

[40] F. Rizzo et al., *The halo bispectrum multipoles in redshift space*, JCAP **01** (2023) 031 [arXiv:2204.13628] [INSPIRE].

[41] M.M. Ivanov et al., *Cosmology with the galaxy bispectrum multipoles: Optimal estimation and application to BOSS data*, Phys. Rev. D **107** (2023) 083515 [arXiv:2302.04414] [INSPIRE].

[42] G. D'Amico et al., *The one-loop bispectrum of galaxies in redshift space from the Effective Field Theory of Large-Scale Structure*, JCAP **07** (2024) 041 [arXiv:2211.17130] [INSPIRE].

[43] L. Verde and A.F. Heavens, *On the trispectrum as a Gaussian test for cosmology*, Astrophys. J. **553** (2001) 14 [astro-ph/0101143] [INSPIRE].

[44] R. Scoccimarro, H.A. Feldman, J.N. Fry and J.A. Frieman, *The bispectrum of IRAS redshift catalogs*, Astrophys. J. **546** (2001) 652 [astro-ph/0004087] [INSPIRE].

[45] H. Gil-Marín et al., *The power spectrum and bispectrum of SDSS DR11 BOSS galaxies — I. Bias and gravity*, Mon. Not. Roy. Astron. Soc. **451** (2015) 539 [arXiv:1407.5668] [INSPIRE].
[46] H. Gil-Marín et al., The power spectrum and bispectrum of SDSS DR11 BOSS galaxies — II. Cosmological interpretation, Mon. Not. Roy. Astron. Soc. 452 (2015) 1914 [arXiv:1408.0027] [INSPIRE].

[47] H. Gil-Marín et al., The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the power spectrum and bispectrum of the DR12 BOSS galaxies, Mon. Not. Roy. Astron. Soc. 465 (2017) 1757 [arXiv:1606.00439] [INSPIRE].

[48] O.H.E. Philcox and M.M. Ivanov, BOSS DR12 full-shape cosmology: $\Lambda$CDM constraints from the large-scale galaxy power spectrum and bispectrum monopole, Phys. Rev. D 105 (2022) 043517 [arXiv:2112.04515] [INSPIRE].

[49] G. D'Amico et al., The BOSS bispectrum analysis at one loop from the Effective Field Theory of Large-Scale Structure, JCAP 05 (2024) 059 [arXiv:2206.08327] [INSPIRE].

[50] G. D'Amico, M. Lewandowski, L. Senatore and P. Zhang, Limits on primordial non-Gaussianities from BOSS galaxy-clustering data, Phys. Rev. D 111 (2025) 063514 [arXiv:2201.11518] [INSPIRE].

[51] D. Gualdi et al., Geometrical compression: a new method to enhance the BOSS galaxy bispectrum monopole constraints, Mon. Not. Roy. Astron. Soc. 484 (2019) L29 [arXiv:1901.00987] [INSPIRE].

[52] D. Gualdi et al., Enhancing BOSS bispectrum cosmological constraints with maximal compression, Mon. Not. Roy. Astron. Soc. 484 (2019) 3713 [arXiv:1806.02853] [INSPIRE].

[53] G. Cabass et al., Constraints on Single-Field Inflation from the BOSS Galaxy Survey, Phys. Rev. Lett. 129 (2022) 021301 [arXiv:2201.07238] [INSPIRE].

[54] G. Cabass et al., Constraints on multifield inflation from the BOSS galaxy survey, Phys. Rev. D 106 (2022) 043506 [arXiv:2204.01781] [INSPIRE].

[55] DESI collaboration, Data Release 1 of the Dark Energy Spectroscopic Instrument, arXiv:2503.14745 [INSPIRE].

[56] S. Novell-Masot, D. Gualdi, H. Gil-Marín and L. Verde, GEO-FPT: a model of the galaxy bispectrum at mildly non-linear scales, JCAP 11 (2023) 044 [arXiv:2303.15510] [INSPIRE].

[57] H. Gil-Marín et al., Perturbation theory approach for the power spectrum: from dark matter in real space to haloes in redshift space, JCAP 11 (2012) 029 [arXiv:1209.3771] [INSPIRE].

[58] BOSS collaboration, The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the LOS-dependent power spectrum of DR12 BOSS galaxies, Mon. Not. Roy. Astron. Soc. 460 (2016) 4188 [arXiv:1509.06386] [INSPIRE].

[59] EBOSS collaboration, The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: measurement of the BAO and growth rate of structure of the luminous red galaxy sample from the anisotropic power spectrum between redshifts 0.6 and 1.0, Mon. Not. Roy. Astron. Soc. 498 (2020) 2492 [arXiv:2007.08994] [INSPIRE].

[60] S. Brieden, H. Gil-Marín and L. Verde, ShapeFit: extracting the power spectrum shape information in galaxy surveys beyond BAO and RSD, JCAP 12 (2021) 054 [arXiv:2106.07641] [INSPIRE].

[61] DESI collaboration, DESI 2024 V: Full-Shape Galaxy Clustering from Galaxies and Quasars, arXiv:2411.12021 [INSPIRE].
[62] EBOSS collaboration, *Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory*, Phys. Rev. D 103 (2021) 083533 [arXiv:2007.08991] [INSPIRE].

[63] C.H. Hahn et al., *The DESI Bright Galaxy Survey: Final Target Selection, Design, and Validation*, Astron. J. 165 (2023) 253 [arXiv:2208.08512] [INSPIRE].

[64] A. Raichoor et al., *Target Selection and Validation of DESI Emission Line Galaxies*, Astron. J. 165 (2023) 126 [arXiv:2208.08513] [INSPIRE].

[65] E. Chaussidon et al., *Target Selection and Validation of DESI Quasars*, Astrophys. J. 944 (2023) 107 [arXiv:2208.08511] [INSPIRE].

[66] DESI collaboration, *DESI 2024 VII: Cosmological Constraints from the Full-Shape Modeling of Clustering Measurements*, arXiv:2411.12022 [INSPIRE].

[67] R. Zhao et al., *Impact and mitigation of imaging systematics for DESI 2024 full shape analysis*, in preparation.

[68] M. Tegmark, *Measuring cosmological parameters with galaxy surveys*, Phys. Rev. Lett. 79 (1997) 3806 [astro-ph/9706198] [INSPIRE].

[69] H.A. Feldman, N. Kaiser and J.A. Peacock, *Power spectrum analysis of three-dimensional redshift surveys*, Astrophys. J. 426 (1994) 23 [astro-ph/9304022] [INSPIRE].

[70] S. Brieden, H. Gil-Marín, L. Verde and J.L. Bernal, *Blind Observers of the Sky*, JCAP 09 (2020) 052 [arXiv:2006.10857] [INSPIRE].

[71] E. Chaussidon et al., *Blinding scheme for the scale-dependence bias signature of local primordial non-Gaussianity for DESI 2024*, JCAP 01 (2025) 135 [arXiv:2406.00191] [INSPIRE].

[72] S. Brieden, H. Gil-Marín and L. Verde, *Model-independent versus model-dependent interpretation of the SDSS-III BOSS power spectrum: Bridging the divide*, Phys. Rev. D 104 (2021) L121301 [arXiv:2106.11931] [INSPIRE].

[73] DESI collaboration, *Validating the galaxy and quasar catalog-level blinding scheme for the DESI 2024 analysis*, JCAP 01 (2025) 128 [arXiv:2404.07282] [INSPIRE].

[74] S. Novell-Masot et al., *Catalog-level blinding on the bispectrum for DESI-like galaxy surveys*, JCAP 10 (2024) 089 [arXiv:2407.12931] [INSPIRE].

[75] N.A. Maksimova et al., *AbacusSummit: a massive set of high-accuracy, high-resolution N-body simulations*, Mon. Not. Roy. Astron. Soc. 508 (2021) 4017 [arXiv:2110.11398] [INSPIRE].

[76] S. Yuan et al., *AbacusHOD: a highly efficient extended multitracer HOD framework and its application to BOSS and eBOSS data*, Mon. Not. Roy. Astron. Soc. 510 (2022) 3301 [arXiv:2110.11412] [INSPIRE].

[77] S. Yuan et al., *The DESI one-per cent survey: exploring the halo occupation distribution of luminous red galaxies and quasi-stellar objects with AbacusSummit*, Mon. Not. Roy. Astron. Soc. 530 (2024) 947 [arXiv:2306.06314] [INSPIRE].

[78] DESI collaboration, *Production of alternate realizations of DESI fiber assignment for unbiased clustering measurement in data and simulations*, JCAP 01 (2025) 127 [arXiv:2404.03006] [INSPIRE].

[79] S. Novell-Masot et al., In preparation.

[80] S. Brieden, H. Gil-Marín and L. Verde, *PT challenge: validation of ShapeFit on large-volume, high-resolution mocks*, JCAP 06 (2022) 005 [arXiv:2201.08400] [INSPIRE].
[81] S. Brieden, H. Gil-Marín and L. Verde, *Model-agnostic interpretation of 10 billion years of cosmic evolution traced by BOSS and eBOSS data*, JCAP 08 (2022) 024 [arXiv:2204.11868] [INSPIRE].

[82] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: baryon acoustic oscillations in the Fourier space*, Mon. Not. Roy. Astron. Soc. 464 (2017) 3409 [arXiv:1607.03149] [INSPIRE].

[83] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: On the measurement of growth rate using galaxy correlation functions*, Mon. Not. Roy. Astron. Soc. 469 (2017) 1369 [arXiv:1607.03148] [INSPIRE].

[84] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: Cosmological implications of the Fourier space wedges of the final sample*, Mon. Not. Roy. Astron. Soc. 467 (2017) 2085 [arXiv:1607.03143] [INSPIRE].

[85] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological implications of the configuration-space clustering wedges*, Mon. Not. Roy. Astron. Soc. 464 (2017) 1640 [arXiv:1607.03147] [INSPIRE].

[86] EBOSS collaboration, *The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: BAO and RSD measurements from anisotropic clustering analysis of the Quasar Sample in configuration space between redshift 0.8 and 2.2*, Mon. Not. Roy. Astron. Soc. 500 (2020) 1201 [arXiv:2007.08998] [INSPIRE].

[87] EBOSS collaboration, *The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: measurement of the BAO and growth rate of structure of the luminous red galaxy sample from the anisotropic correlation function between redshifts 0.6 and 1*, Mon. Not. Roy. Astron. Soc. 500 (2020) 736 [arXiv:2007.08993] [INSPIRE].

[88] EBOSS collaboration, *The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: BAO and RSD measurements from the anisotropic power spectrum of the quasar sample between redshift 0.8 and 2.2*, Mon. Not. Roy. Astron. Soc. 499 (2020) 210 [arXiv:2007.08999] [INSPIRE].

[89] EBOSS collaboration, *The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: measurement of the BAO and growth rate of structure of the emission line galaxy sample from the anisotropic power spectrum between redshift 0.6 and 1.1*, Mon. Not. Roy. Astron. Soc. 501 (2021) 5616 [arXiv:2007.09008] [INSPIRE].

[90] C. Alcock and B. Paczynski, *An evolution free test for non-zero cosmological constant*, Nature 281 (1979) 358 [INSPIRE].

[91] H.E. Noriega et al., *Comparing Compressed and Full-Modeling analyses with FOLPS: implications for DESI 2024 and beyond*, JCAP 01 (2025) 136 [arXiv:2404.07269] [INSPIRE].

[92] T. Baldauf, U. Seljak, V. Desjacques and P. McDonald, *Evidence for Quadratic Tidal Tensor Bias from the Halo Bispectrum*, Phys. Rev. D 86 (2012) 083540 [arXiv:1201.4827] [INSPIRE].

[93] S. Saito et al., *Understanding higher-order nonlocal halo bias at large scales by combining the power spectrum with the bispectrum*, Phys. Rev. D 90 (2014) 123522 [arXiv:1405.1447] [INSPIRE].

[94] M. Crocce and R. Scoccimarro, *Renormalized cosmological perturbation theory*, Phys. Rev. D 73 (2006) 063519 [astro-ph/0509418] [INSPIRE].

[95] M. Maus et al., *A comparison of effective field theory models of redshift space galaxy power spectra for DESI 2024 and future surveys*, JCAP 01 (2025) 134 [arXiv:2404.07272] [INSPIRE].
[96] M. Maus et al., *An analysis of parameter compression and Full-Modeling techniques with Velocileptors for DESI 2024 and beyond*, JCAP 01 (2025) 138 [arXiv:2404.07312] [INSPIRE].

[97] Y. Lai et al., *A comparison between ShapeFit compression and Full-Modelling method with PyBird for DESI 2024 and beyond*, JCAP 01 (2025) 139 [arXiv:2404.07283] [INSPIRE].

[98] DESI collaboration, *Full Modeling and parameter compression methods in configuration space for DESI 2024 and beyond*, JCAP 01 (2025) 129 [arXiv:2404.07268] [INSPIRE].

[99] M.J. Wilson, J.A. Peacock, A.N. Taylor and S. de la Torre, *Rapid modelling of the redshift-space power spectrum multipoles for a masked density field*, Mon. Not. Roy. Astron. Soc. 464 (2017) 3121 [arXiv:1511.07799] [INSPIRE].

[100] M. Pinon et al., *Mitigation of DESI fiber assignment incompleteness effect on two-point clustering with small angular scale truncated estimators*, JCAP 01 (2025) 131 [arXiv:2406.04804] [INSPIRE].

[101] M. Maus, S.-F. Chen and M. White, *A comparison of template vs. direct model fitting for redshift-space distortions in BOSS*, JCAP 06 (2023) 005 [arXiv:2302.07430] [INSPIRE].

[102] R. Scoccimarro, H.M.P. Couchman and J.A. Frieman, *The Bispectrum as a Signature of Gravitational Instability in Redshift-Space*, Astrophys. J. 517 (1999) 531 [astro-ph/9808305] [INSPIRE].

[103] N.S. Sugiyama, S. Saito, F. Beutler and H.-J. Seo, *A complete FFT-based decomposition formalism for the redshift-space bispectrum*, Mon. Not. Roy. Astron. Soc. 484 (2019) 364 [arXiv:1803.02132] [INSPIRE].

[104] M.S. Wang et al., *Window convolution of the galaxy clustering bispectrum*, arXiv:2411.14947 [INSPIRE].

[105] K. Pardede et al., *Bispectrum-window convolution via Hankel transform*, JCAP 10 (2022) 066 [arXiv:2203.04174] [INSPIRE].

[106] N. Metropolis et al., *Equation of state calculations by fast computing machines*, J. Chem. Phys. 21 (1953) 1087 [INSPIRE].

[107] W.K. Hastings, *Monte Carlo Sampling Methods Using Markov Chains and Their Applications*, Biometrika 57 (1970) 97 [INSPIRE].

[108] S. Novell-Masot, H. Gil-Marín and L. Verde, *On approximations of the redshift-space bispectrum and power spectrum multipoles covariance matrix*, JCAP 06 (2024) 048 [arXiv:2306.03137] [INSPIRE].

[109] J. Hartlap, P. Simon and P. Schneider, *Why your model parameter confidences might be too optimistic: Unbiased estimation of the inverse covariance matrix*, Astron. Astrophys. 464 (2007) 399 [astro-ph/0608064] [INSPIRE].

[110] E. Sellentin and A.F. Heavens, *Parameter inference with estimated covariance matrices*, Mon. Not. Roy. Astron. Soc. 456 (2016) L132 [arXiv:1511.05969] [INSPIRE].

[111] A.J. Ross et al., *The construction of large-scale structure catalogs for the Dark Energy Spectroscopic Instrument*, JCAP 01 (2025) 125 [arXiv:2405.16593] [INSPIRE].

[112] N. Findlay et al., *Exploring HOD-dependent systematics for the DESI 2024 Full-Shape galaxy clustering analysis*, arXiv:2411.12023 [INSPIRE].

[113] DESI collaboration, *HOD-dependent systematics for luminous red galaxies in the DESI 2024 BAO analysis*, JCAP 01 (2025) 133 [arXiv:2404.03008] [INSPIRE].
[114] R. Gsponer et al., *Fiducial cosmology impact for DESI 2024 full shape analysis*, in preparation.

[115] S. Sanz-Wuhl, H. Gil-Marín, A.J. Cuesta and L. Verde, *BAO cosmology in non-spatially flat background geometry from BOSS+eBOSS and lessons for future surveys*, JCAP 05 (2024) 116 [arXiv:2402.03427] [INSPIRE].

[116] DESI collaboration, *Fiducial-cosmology-dependent systematics for the DESI 2024 BAO analysis*, JCAP 01 (2025) 144 [arXiv:2406.06085] [INSPIRE].

[117] E.V. Linder, *Cosmic growth history and expansion history*, Phys. Rev. D 72 (2005) 043529 [astro-ph/0507263] [INSPIRE].

[118] E.V. Linder and R.N. Cahn, *Parameterized Beyond-Einstein Growth*, Astropart. Phys. 28 (2007) 481 [astro-ph/0701317] [INSPIRE].

[119] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, *emcee: The MCMC Hammer*, Publ. Astron. Soc. Pac. 125 (2013) 306 [arXiv:1202.3665] [INSPIRE].

[120] M.C. Galassi et al., *GNU Scientific Library*, Network Theory, Ltd. (2019) [INSPIRE].

[121] P. Virtanen et al., *SciPy 1.0 — Fundamental Algorithms for Scientific Computing in Python*, Nature Meth. 17 (2020) 261 [arXiv:1907.10121] [INSPIRE].

[122] C.R. Harris et al., *Array programming with NumPy*, Nature 585 (2020) 357 [arXiv:2006.10256] [INSPIRE].

[123] A. Lewis, *GetDist: a Python package for analysing Monte Carlo samples*, arXiv:1910.13970 [INSPIRE].

[124] ASTROPY collaboration, *The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package, Astrophys. J.* 935 (2022) 167 [arXiv:2206.14220] [INSPIRE].

[125] J.D. Hunter, *Matplotlib: A 2D Graphics Environment*, Comput. Sci. Eng. 9 (2007) 90 [INSPIRE].

[126] H. Gil-Marín et al., *An improved fitting formula for the dark matter bispectrum*, JCAP 02 (2012) 047 [arXiv:1111.4477] [INSPIRE].

[127] A. Taruya, T. Nishimichi and S. Saito, *Baryon Acoustic Oscillations in 2D: Modeling Redshift-space Power Spectrum from Perturbation Theory*, Phys. Rev. D 82 (2010) 063522 [arXiv:1006.0699] [INSPIRE].

[128] T. Nishimichi and A. Taruya, *Baryon Acoustic Oscillations in 2D II: Redshift-space halo clustering in N-body simulations*, Phys. Rev. D 84 (2011) 043526 [arXiv:1106.4562] [INSPIRE].

[129] BOSS collaboration, *The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Testing gravity with redshift-space distortions using the power spectrum multipoles*, Mon. Not. Roy. Astron. Soc. 443 (2014) 1065 [arXiv:1312.4611] [INSPIRE].

[130] K.C. Chan, R. Scoccimarro and R.K. Sheth, *Gravity and Large-Scale Non-local Bias*, Phys. Rev. D 85 (2012) 083509 [arXiv:1201.3614] [INSPIRE].

[131] J.C. Jackson, *A critique of Rees’s theory of primordial gravitational radiation*, Mon. Not. Roy. Astron. Soc. 156 (1972) 1P [arXiv:0810.3908] [INSPIRE].

[132] J. Fry, *Gravity, bias, and the galaxy three-point correlation function*, Phys. Rev. Lett. 73 (1994) 215.

[133] R. Scoccimarro et al., *Nonlinear evolution of the bispectrum of cosmological perturbations*, Astrophys. J. 496 (1998) 586 [astro-ph/9704075] [INSPIRE].
[134] H. Gil-Marín et al., Dark matter and halo bispectrum in redshift space: theory and applications, JCAP 12 (2014) 029 [arXiv:1407.1836] [INSPIRE].

[135] A.J.S. Hamilton, Measuring Omega and the real correlation function from the redshift correlation function, Astrophys. J. Lett. 385 (1992) L5 [INSPIRE].

[136] S. Cole, K.B. Fisher and D.H. Weinberg, Fourier analysis of redshift-space distortions and the determination of $\omega$, Mon. Not. Roy. Astron. Soc. 267 (1994) 785 [astro-ph/9308003] [INSPIRE].
Author List

S. Novell-Masot ⋆ 1, H. Gil-Marín ⋆ 2,3,1, L. Verde ⋆ 4,1, J. Aguilar 5, S. Ahlen ⋆ 6, S. Bailey ⋆ 5,

S. BenZvi ⋆ 7, D. Bianchi ⋆ 8,9, D. Brooks 10, E. Buckley-Geer 11,12, A. Carnero Rosell ⋆ 13,14,

E. Chaussidon ⋆ 5, T. Claybaugh 5, S. Cole ⋆ 15, A. Cuceu ⋆ 5, K.S. Dawson ⋆ 16, A. de la Macorra ⋆ 17,

R. Demina 7, A. Dey ⋆ 18, B. Dey ⋆ 19,20, P. Doel 10, S. Ferraro ⋆ 5,21, A. Font-Ribera ⋆ 22,

J.E. Forero-Romero ⋆ 23,24, E. Gaztañaga 3,25,26, S. Gontcho A Gontcho ⋆ 5,

A.X. Gonzalez-Morales ⋆ 27, G. Gutierrez 12, H.K. Herrera-Alcantar ⋆ 28,29, K. Honscheid ⋆ 30,31,32,

C. Howlett ⋆ 33, S. Juneau ⋆ 18, R. Kehoe 34, D. Kirkby ⋆ 35, T. Kisner ⋆ 5, A. Kremin ⋆ 5,

C. Lamman ⋆ 36, M. Landriau ⋆ 5, L. Le Guillou ⋆ 37, M.E. Levi ⋆ 5, C. Magneville 29,

M. Manera ⋆ 38,22, A. Meisner ⋆ 18, R. Miquel 4,22, J. Moustakas ⋆ 39, A. Muñoz-Gutiérrez 17,

A.D. Myers 40, S. Nadathur ⋆ 25, G. Niz ⋆ 27,41, H.E. Noriega ⋆ 42,17, W.J. Percival ⋆ 43,44,45,

C. Poppett 5,46,21, F. Prada ⋆ 47, I. Pérez-Ràfols ⋆ 48, A.J. Ross ⋆ 30,49,32, G. Rossi 50,

L. Samushia ⋆ 51,52,53, E. Sanchez ⋆ 54, D. Schlegel 5, M. Schubnell 55,56, H. Seo ⋆ 57, J. Silber ⋆ 5,

D. Sprayberry 18, G. Tarlé ⋆ 56, M. Vargas-Magaña ⋆ 17, B.A. Weaver 18, P. Zarrouk ⋆ 37, R. Zhou ⋆ 5,

H. Zou ⋆ 58

1 Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB),
c. Martí i Franquès, 1, 08028 Barcelona, Spain

2 Departament de Física Quàntica i Astrofísica, Universitat de Barcelona,
Martí i Franquès 1, E08028 Barcelona, Spain

3 Institut d'Estudis Espacials de Catalunya (IEEC),
c/ Esteve Terradas 1, Edifici RDIT, Campus PMT-UPC, 08860 Castelldefels, Spain

4 Institució Catalana de Recerca i Estudis Avançats,
Passeig de Lluís Companys, 23, 08010 Barcelona, Spain

5 Lawrence Berkeley National Laboratory,
1 Cyclotron Road, Berkeley, CA 94720, U.S.A.

6 Physics Dept., Boston University,
590 Commonwealth Avenue, Boston, MA 02215, U.S.A.

7 Department of Physics & Astronomy, University of Rochester,
206 Bausch and Lomb Hall, P.O. Box 270171, Rochester, NY 14627-0171, U.S.A.

8 Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano,
Via Celoria 16, I-20133 Milano, Italy

9 INAF-Osservatorio Astronomico di Brera,
Via Brera 28, 20122 Milano, Italy

10 Department of Physics & Astronomy, University College London,
Gower Street, London, WC1E 6BT, U.K.

11 Department of Astronomy and Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, IL 60637, U.S.A.

12 Fermi National Accelerator Laboratory,
PO Box 500, Batavia, IL 60510, U.S.A.

13 Departamento de Astrofísica, Universidad de La Laguna (ULL),
E-38206, La Laguna, Tenerife, Spain

14 Instituto de Astrofísica de Canarias,
C/ Vía Láctea, s/n, E-38205 La Laguna, Tenerife, Spain

15 Institute for Computational Cosmology, Department of Physics, Durham University,
South Road, Durham DH1 3LE, U.K.

16 Department of Physics and Astronomy, The University of Utah,
115 South 1400 East, Salt Lake City, UT 84112, U.S.A.

17 Instituto de Física, Universidad Nacional Autónoma de México,
Circuito de la Investigación Científica, Ciudad Universitaria,
Cd. de México C.P. 04510, México
¹⁸ NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, U.S.A.

¹⁹ Department of Astronomy & Astrophysics, University of Toronto,
Toronto, ON M5S 3H4, Canada

²⁰ Department of Physics & Astronomy and Pittsburgh Particle Physics,
Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh,
3941 O'Hara Street, Pittsburgh, PA 15260, U.S.A.

²¹ University of California,
Berkeley, 110 Sproul Hall #5800 Berkeley, CA 94720, U.S.A.

²² Institut de Física d'Altes Energies (IFAE),
The Barcelona Institute of Science and Technology,
Edifici Cn, Campus UAB, 08193, Bellaterra (Barcelona), Spain

²³ Departamento de Física, Universidad de los Andes,
Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia

²⁴ Observatorio Astronómico, Universidad de los Andes,
Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia

²⁵ Institute of Cosmology and Gravitation, University of Portsmouth,
Dennis Sciama Building, Portsmouth, PO1 3FX, U.K.

²⁶ Institute of Space Sciences, ICE-CSIC,
Campus UAB, Carrer de Can Magrans s/n, 08913 Bellaterra, Barcelona, Spain

²⁷ Departamento de Física, DCI-Campus León, Universidad de Guanajuato,
Loma del Bosque 103, León, Guanajuato C. P. 37150, México

²⁸ Institut d'Astrophysique de Paris. 98 bis boulevard Arago. 75014 Paris, France

²⁹ IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

³⁰ Center for Cosmology and AstroParticle Physics, The Ohio State University,
191 West Woodruff Avenue, Columbus, OH 43210, U.S.A.

³¹ Department of Physics, The Ohio State University,
191 West Woodruff Avenue, Columbus, OH 43210, U.S.A.

³² The Ohio State University, Columbus, 43210 OH, U.S.A.

³³ School of Mathematics and Physics, University of Queensland,
Brisbane, QLD 4072, Australia

³⁴ Department of Physics, Southern Methodist University,
3215 Daniel Avenue, Dallas, TX 75275, U.S.A.

³⁵ Department of Physics and Astronomy, University of California,
Irvine, 92697, U.S.A.

³⁶ Center for Astrophysics | Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, U.S.A.

³⁷ Sorbonne Université, CNRS/IN2P3,
Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE),
FR-75005 Paris, France

³⁸ Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona,
08193 Bellaterra (Barcelona), Spain

³⁹ Department of Physics and Astronomy, Siena College,
515 Loudon Road, Loudonville, NY 12211, U.S.A.

⁴⁰ Department of Physics & Astronomy, University of Wyoming,
1000 E. University, Dept. 3905, Laramie, WY 82071, U.S.A.

⁴¹ Instituto Avanzado de Cosmología A.C.,
San Marcos 11 — Atenas 202, Magdalena Contreras, Ciudad de México C. P. 10720, México

⁴² Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México,
Av. Universidad s/n, Cuernavaca, Morelos, C. P. 62210, México

⁴³ Department of Physics and Astronomy, University of Waterloo,
200 University Ave W, Waterloo, ON N2L 3G1, Canada

⁴⁴ Perimeter Institute for Theoretical Physics,
31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada
<sup>45</sup> Waterloo Centre for Astrophysics, University of Waterloo,
200 University Ave W, Waterloo, ON N2L 3G1, Canada

<sup>46</sup> Space Sciences Laboratory, University of California,
Berkeley, 7 Gauss Way, Berkeley, CA 94720, U.S.A.

<sup>47</sup> Instituto de Astrofísica de Andalucía (CSIC),
Glorieta de la Astronomía, s/n, E-18008 Granada, Spain

<sup>48</sup> Departament de Física, EEBE, Universitat Politècnica de Catalunya,
c/Eduard Maristany 10, 08930 Barcelona, Spain

<sup>49</sup> Department of Astronomy, The Ohio State University,
4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, U.S.A.

<sup>50</sup> Department of Physics and Astronomy, Sejong University,
209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Republic of Korea

<sup>51</sup> Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia

<sup>52</sup> Department of Physics, Kansas State University,
116 Cardwell Hall, Manhattan, KS 66506, U.S.A.

<sup>53</sup> Faculty of Natural Sciences and Medicine, Ilia State University,
0194 Tbilisi, Georgia

<sup>54</sup> CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

<sup>55</sup> Department of Physics, University of Michigan,
450 Church Street, Ann Arbor, MI 48109, U.S.A.

<sup>56</sup> University of Michigan, 500 S. State Street, Ann Arbor, MI 48109, U.S.A.

<sup>57</sup> Department of Physics & Astronomy, Ohio University,
139 University Terrace, Athens, OH 45701, U.S.A.

<sup>58</sup> National Astronomical Observatories, Chinese Academy of Sciences,
A20 Datun Rd., Chaoyang District, Beijing, 100012, P.R. China