DESI DR2 Results II: Measurements of Baryon Acoustic Oscillations and Cosmological Constraints

(DESI Collaboration)* {sup}21` Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA

{sup}22` CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182 21, Prague, Czech Republic.

{sup}23` Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, UK

{sup}24` Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Edifici Cn, Campus UAB, 08193, Bellaterra (Barcelona), Spain

{sup}25` Institut d’Estudis Espacials de Catalunya (IEEC), c/ Esteve Terradas 1, Edifici RDIT, Campus PMT-UPC, 08860 Castelldefels, Spain

{sup}26` Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans s/n, 08913 Bellaterra, Barcelona, Spain

{sup}27` Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, U.K

{sup}28` The Ohio State University, Columbus, 43210 OH, USA

{sup}29` Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK

{sup}30` Institute of Astronomy and Department of Physics, National Tsing Hua University, 101 Kuang-Fu Rd. Sec. 2, Hsinchu 30013, Taiwan

{sup}31` NASA Einstein Fellow

{sup}32` Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

{sup}33` 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

{sup}34` CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

{sup}35` Department of Astronomy, San Diego State University, 5500 Campanile Drive, San Diego, CA 92182, USA

{sup}36` Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

{sup}37` Department of Physics & Astronomy and Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O’Hara Street, Pittsburgh, PA 15260, USA

{sup}38` University of Chinese Academy of Sciences, Nanjing 211135, People’s Republic of China.

{sup}39` Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA

{sup}40` University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA 94720, USA

{sup}41` Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

{sup}42` Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Menlo Park, CA 94305, USA

{sup}43` SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA

{sup}44` Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain

{sup}45` Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio IP, CP 111711, Bogotá, Colombia

{sup}46` Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio IH, CP 111711 Bogotá, Colombia

{sup}47 Center for Computational Astrophysics, Flatiron Institute, 162 5{sup th} Avenue, New York, NY 10010, USA

{sup}48 Scientific Computing Core, Flatiron Institute, 162 5{sup th} Avenue, New York, NY 10010, USA

{sup}49` Department de Física Quàntica i Astrofísica, Universitat de Barcelona, Martí i Franquès 1, E08028 Barcelona, Spain

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

{sup}51` Departamento de Física, DCI-Campus León, Universidad de Guanajuato, Loma del Bosque 103, León, Guanajuato C. P. 37150, México

{sup}52` Department of Physics and Astronomy, University of California, Irvine, 92697, USA

{sup}53` Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

{sup}54` Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA

{sup}55` Institute of Physics, Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, Chemin Pegasi 51, CH-1290 Versoix, Switzerland

{sup}56` Institut d’Astrophysique de Paris, 98 bis boulevard Arago, 75014 Paris, France

{sup}57` Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA

{sup}58` Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA

{sup}59` School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

{sup}60` Department of Physics, The University of Texas at Dallas, 800 W. Campbell Rd., Richardson, TX 75080, USA

{sup}61` Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA

{sup}62` Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA

{sup}63` Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada

{sup}64` Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada

{sup}65` Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada

{sup}66` Astrophysics & Space Institute, Schmidt Sciences, New York, NY 10011, USA

{sup}67` Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), FR-75005 Paris, France

{sup}68` Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California, 1156 High Street, Santa Cruz, CA 95064, USA

{sup}69` Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95065, USA

{sup}70` Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea

{sup}71` University of Science and Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon 34113, Republic of Korea

{sup}72` Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Spain

• [email protected]

Introduction

Cosmic acceleration remains the most pressing problem in contemporary cosmology, implying a pervasive new form of energy with exotic physical properties, or a breakdown of Einstein gravity on cosmological scales, or perhaps both. To probe the physics of acceleration, cosmologists seek to measure the history of cosmic expansion and the history of gravitational clustering with the greatest achievable precision over a wide span of redshift. Baryon acoustic oscillations (BAO) provide a powerful tool for measuring the expansion history, using a characteristic scale that is imprinted on matter clustering by pressure waves that propagate in the coupled baryon-photon fluid of the pre-recombination Universe [1, 2, 3]. This paper examines the cosmological implications of the BAO measurements from the second data release (DR2) of the Dark Energy Spectroscopic Instrument (DESI, [4, 5, 6]), consisting of data from the first three years of operation.

Since the first clear detections of BAO in the Sloan Digital Sky Survey and the Two-Degree Field Galaxy Redshift Survey [7, 8], BAO measurements have played a central role in observational cosmology. The key technical requirement is a large-volume spectroscopic survey with sufficient sampling density, and previous surveys designed with BAO measurements as a defining goal include WiggleZ [9], the Baryon Oscillation Spectroscopic Survey (BOSS) [10] of SDSS-III [11], and its extension eBOSS [12] in SDSS-IV [13]. In addition to galaxy and quasar redshifts, BOSS and eBOSS measured BAO in the Ly$\alpha$ forest absorption spectra of $z > 2$ quasars, an approach first proposed by [14, 15]. Transverse BAO can also be measured in photometric surveys (e.g., [16]), though the precision obtained for a given number of tracers is much higher with spectroscopic redshifts. DESI is designed specifically to enable a spectroscopic BAO survey of unprecedented power and efficiency [17], as shown in its survey validation [18] based on the early data release [19]. In its first year of observations [20], DESI already achieved BAO measurements competitive with those of all previous surveys combined [21, 22]. Now, with redshifts of more than 30 million galaxies and quasars, and Ly$\alpha$ forest spectra of more than 820,000 quasars, the DESI DR2 sample is by far the largest spectroscopic galaxy sample to date.

The BAO technique (reviewed in §4 of [23]), which is key to DESI, exploits the enhancement of clustering at the scale of the pre-recombination sound horizon,

Here $c_s(z)$ is the speed of sound in the photon-baryon fluid, and $z_d \approx1060$ [24] is the redshift at which acoustic waves stall because photons no longer ‘drag’ the baryons. Assuming standard pre-recombination physics, the sound horizon can be computed given the densities of baryons, cold dark matter (CDM), photons, and other relativistic species [25],

Eq. (2) is scaled to the best-fit values from Planck [24] of $\omega_b \equiv\Omega_b h^2$ and $\omega_{bc} \equiv(\Omega_b+\Omega_c)h^2$ and to the energy content of three neutrino species that are fully relativistic at $z > z_d$. Here $h \equiv H_0/(100\,\mathrm{km\,s^{-1}\,Mpc^{-1}})$ and $\Omega_b$ and $\Omega_c$ denote the present day fractional energy densities relative to critical in baryons and cold dark matter, respectively. A measurement of the BAO scale in the transverse direction at redshift $z$ constrains the transverse comoving distance, which is given by

where $\Omega_K$ is the curvature density parameter, which converges to the flat universe case

in the limit $\lvert\Omega_K\rvert\ll1$. A measurement in the line-of-sight direction constrains the expansion rate $H(z)$ or the corresponding distance,

Because the inferred distances are relative to the sound horizon, the directly constrained quantities are the ratios $D_M/r_d$ and $D_H/r_d$.

The combination of BAO and cosmic microwave background (CMB) anisotropies is powerful for two reasons. First, the CMB provides tight constraints on $\omega_b$ and $\omega_{bc}$, leading to a 0.2% determination of $r_d$ from eq. (2) (for standard $N_{\mathrm{eff}} = 3.04$). As a result, the BAO+CMB combination allows absolute measurements of $D_M(z)$ and $D_H(z)$. Second, the same physics imprints both the BAO and the acoustic peaks in the CMB power spectrum [26]. The angular scale of these peaks, denoted $\theta_*$, is measured with exquisite precision (fractional error $\sim10^{-4}$), giving a near-perfect measurement of $D_M(z_*)/r_*$, where $r_*$ is the comoving sound horizon at the end of recombination, at redshift $z_* \approx1089$. Figure 1 shows a pedagogical view of the cosmological role of these measurements of the expansion history.

Type Ia supernovae (SNe) are standardizable candles which serve as low-redshift probes in addition to BAO, measuring the luminosity distance $D_L(z) = (1+z)D_M(z)$. The SNe Hubble diagram provided the first direct evidence for cosmic acceleration [27, 28], and subsequent cosmological surveys have obtained well measured light curves for many hundreds of supernovae out to $z = 1$ and beyond (e.g., [29, 30, 31], and references therein). BAO and SNe measurements constrain dark energy and neutrino masses through their impact on the background evolution, thus determining $H(z)$. Assuming that general relativity (GR) correctly describes the dynamics of expansion, the evolution of $H(z)$ is governed by the Friedmann equation, which can be written

Here $\Omega_{bc} = \Omega_b + \Omega_c$ and $\Omega_\gamma$, $\Omega_K$, $\Omega_\nu$ and $\Omega_{\mathrm{DE}}$ refer to the energy densities in radiation, curvature, neutrinos and dark energy, respectively. We refer to the fractional energy density in matter as $\Omega_m$, which includes neutrinos when they are non-relativistic. The energy densities of baryons and cold dark matter scale as $(1+z)^3$. The scaling of the neutrino energy density transitions from $(1+z)^3$ to $(1+z)^4$ at high redshifts, when $(1+z) \sim(m_\nu/5 \times10^{-4}\,\mathrm{eV})$ [32, 33]. The sum of neutrino masses determines the present day density [34]

For BAO cosmology, the important characteristic of ‘CDM’ is that its energy density scales as $(1+z)^3$, both before and after recombination, and that it does not couple non-gravitationally to photons or baryons in a way that affects the scale of the acoustic oscillations. Some variations such as self-interaction or ‘warm’ thermal velocities would affect small scale clustering and galaxy formation but not BAO. Decaying dark matter, on the other hand, would alter the $(1+z)^3$ energy scaling even if the dark matter is cold and non-interacting, thus affecting the BAO scale.

The $\Lambda$CDM model assumes a cosmological constant dark energy ($\Lambda$) with energy density $c^2\rho_{\mathrm{DE}}$ that is constant in space and time. If dark energy has an equation-of-state parameter $w(z) \equiv P(z)/(c^2\rho_{\mathrm{DE}}(z))$ where $P(z)$ is its pressure, then its energy density evolves as

For constant $w$, the r.h.s. of eq. (8) is simply $(1+z)^{3(1+w)}$, and a cosmological constant corresponds to $w=-1$. A commonly used parametric model expresses $w(z)$ in terms of the expansion factor $a=(1+z)^{-1}$,

so that $w$ evolves from a value $\sim(w_0+w_a)$ at high redshift to a present-day value of $w_0$. This parametrization accurately represents the behavior of many physically motivated dark energy models [35, 36, 37], though more complicated evolution is possible. We refer to models that assume CDM and eq. (9) as $w_0w_a$CDM and models with constant $w$ (i.e., $w_a=0$) as $w$CDM. In the $w_0w_a$CDM model the integral of eq. (8) can be evaluated analytically, yielding

Through most of this paper we will assume a flat universe and thus $\Omega_K=0$ in eq. (6), motivated by the tight constraints obtained on $\Omega_K$ when it is allowed to vary freely [38].

The combination of BAO, CMB, and SNe data has allowed tight constraints on the energy density $\Omega_{\mathrm{DE}}$ and equation-of-state $w(z)$ of dark energy, on space curvature $\Omega_K$, on neutrino masses $\sum m_\nu$, and on many possible departures from standard cosmology (see, e.g., [39, 40, 41, 42]). For reviews that explain the complementary constraining power of CMB, BAO, SNe, weak lensing, and other cosmological measurements, see [23, 43, 44]. The analysis of the DESI DR1 measurements from [21, 22] in [38] provided tight constraints on parameters of the $\Lambda$CDM model and intriguing hints of evolving dark energy, with significances ranging from $\sim2.5\sigma$ to $\sim3.9\sigma$ depending on the combination of datasets used for the analysis. The analysis of the full shape of the power spectrum measured with galaxies and quasars [45] [46] confirmed these findings and added new information on the amplitude of perturbations and modified gravity parameters [47].

Rapidly evolving dark energy, with $|w_a| \sim1$, would be an astounding discovery, and these results have inspired both enthusiastic theorizing and healthy skepticism. In our analysis here of the DESI DR2 BAO results, we pay particular attention to the nature and statistical significance of the evidence for evolving dark energy and to how that evidence depends on the choice of datasets. We also examine the constraints on $\sum m_\nu$ from the DESI DR2 data in combination with CMB and SNe, for both $\Lambda$CDM and $w_0w_a$CDM. When we refer to ‘DESI’ alone in tables and figure legends, we treat the BAO as an uncalibrated standard ruler. In some of our $\Lambda$CDM analyses, we examine constraints that adopt a big bang nucleosynthesis (BBN) prior on $\omega_b$, with the value of $\omega_{bc}$ in eq. (2) coming from the model fit itself. We achieve tighter constraints and sharper tests by combining DESI with CMB data that directly constrain $\omega_b$ and $\omega_{bc}$ and add the precise measurement of $\theta_\ast$ at $z = z_\ast$.

This work is accompanied by a set of supporting papers, highlighted in Table I. The structure of this paper is as follows. In § 2, we describe the DESI DR2 data and large-scale structure catalogs. § 3 presents the DESI DR2 distance measurements and internal consistency checks, and presents a comparison with SDSS. In Section IV, we describe the external datasets that will be combined with DESI BAO. § 4.3 introduces our cosmological inference method. Section VI presents cos-

caption: Supporting papers relevant to this work and the corresponding sections where their results are discussed.

mological parameter constraints in a $\Lambda$CDM framework. Section VII presents constraints in a more generalized dark energy framework and analyzes tensions with respect to the $\Lambda$CDM model. Section VIII presents constraints on the neutrino sector and explores the role of neutrino masses in our results. Finally, Section IX presents a summary and our main conclusions.

DESI Data

The DESI Collaboration has measured redshifts for over 30 million galaxies and quasars in just three years of operation, $\sim$14 million of which are used¹ in this analysis, as described below. This extraordinarily high rate of data collection is possible because we built an extremely efficient instrument that can measure thousands of spectra in a single observation [6], combined with the light-gathering power of the 4-m Nicholas U. Mayall Telescope at the Kitt Peak National Observatory. DESI collects the light for 5000 spectra per observation with a robotic focal plane assembly [51] that can quickly align the positions of fiber optics cables [52] across the seven square degree field of view of the prime focus corrector [53]. For each observation, there is a custom focal plane configuration or ‘tile’ that defines the DESI ‘target’ [54] associated with each robotic positioner. When repeated observations of a tile are required to obtain the minimum effective observing time [55], these maintain the same configuration. The light of each of these targets, along with calibration stars and sky spectra, are recorded from 360–980 nm with ten bench-mounted spectrographs that are located in a climate-controlled enclosure. This configuration helps to enable the superb wavelength and flux calibration of the DESI data.

The DESI main survey started observations on 14 May 2021 after a period of survey validation [18]. This paper presents the analysis of the main survey data that will be released with the second data release or DR2,

which includes observations through 9 April 2024. The DESI spectroscopic reduction [56] and redshift estimation (Redrock [57, 58]) pipelines were applied to the DR2 dataset in a homogeneous processing run denoted as ‘Kibo’. An error in the processing involved in the co-addition of spectra from separate exposures was subsequently identified and fixed, and the pipeline was rerun and denoted ‘Loa’. Approximately 0.1% of the measured redshifts change significantly between Kibo and Loa. The LSS catalogs for DR2 BAO measurements were produced for both Kibo and Loa. Decisions that affect the masking of the data (described further below) were made using the Kibo data and were not reconsidered with Loa. Both datasets will be released with DR2.

The DESI survey has two main observing programs, ‘bright’ and ‘dark,’ that are defined based on the nighttime sky conditions, and each of these programs has distinct target classes [54, 55, 59, 60, 61, 62]. The extragalactic sample for the bright program is the ‘bright galaxy sample’ (BGS).² Luminous red galaxies (LRGs), emission line galaxies (ELGs), and quasars (QSOs) are all observed during dark conditions. Dark- and bright-time targets are processed separately. At $z < 2.1$, we measure BAO with the autocorrelation of the confirmed members of each target class, processed through ‘large-scale structure’ (LSS) catalogs as described in § 3 A. At higher redshifts, we measure the auto-correlation of the Ly$\alpha$ forest absorption in the spectra of quasars and the cross-correlation of the forest absorption with quasar positions. The data samples used for these Ly$\alpha$ measurements are described further in § 3 B. The analysis of the Ly$\alpha$ data and BAO results are presented in the companion key paper [63]. DR2 contains 6671 dark and 5171 bright tiles. These are respectively 2.4 times and 2.3 times the number released in DR1.

Galaxy and quasar large-scale Structure Catalogs

For clustering science analyses, the catalogs of measured redshifts and details of the target selection are converted into large-scale structure (LSS) catalogs, as described in [64]. For the specific choices involved in the analysis, we mostly match those decided in [65]. New developments include the following:

• The bad fiber list has been updated, applying the same methods as for DR1 [66], but using the Kibo redshift results.

• For correcting QSO imaging systematics, we switch from a random forest method to the linear regression method used for the LRG and BGS samples, motivated by the conclusions of [67].

• The BGS sample is more dense, as we choose to apply a less restrictive luminosity threshold. We discuss this further below.

Table II contains basic information about the content of the LSS catalogs used for this analysis. In addition to the increase in the raw number of redshifts, one can observe that the sky coverage and the completeness within that sky coverage have increased substantially compared to the DR1 LSS catalogs. Compared to DR1, the size of the ELG sample has grown by a factor of 2.7, while the LRG and QSO samples have grown by factors of 2.1 and 1.7, respectively. The difference in these factors is due to the fact that the coverage, in terms of the typical number of overlapping tiles, has increased substantially between DR1 and DR2. This can be seen by comparing the top panels of Figure 2 to the top panels of the equivalent DR1 figure (figure 2 in [65]). The increase in coverage grows the size of the ELG sample the most, as ELG targets are assigned at the lowest priority (see [54, 64]). Their assignment completeness has increased by 53% (from 35.2% to 53.7%), while the increases are more modest for LRG and QSO (19% and 7%). The increase in area from DR1 is between 54% (QSO) and 75% (LRG and ELG). The areas differ due to the tiles covered (bright vs. dark programs), the priority masking (primarily applied to LRG and ELG due to the influence of QSO) and imaging veto masks (these all differ slightly per tracer). The masks applied match the criteria defined in DR1 [65]. The amount of area in the priority mask is considerably lower than in DR1, despite the greater overall footprint. As was the case for the assignment completeness, this is due to the increased typical number of overlapping tiles: the priority mask is primarily caused by QSO targets receiving higher priority, which happens in the first observation of a tile but not in subsequent overlapping tiles except for quasars at $z > 2.1$. The result is that the area in the LRG and ELG priority mask has reduced from 1667 to 1153 deg$^{2}$, thus narrowing the gap between the size of the QSO and LRG/ELG footprints. The remaining differences between the footprints of the dark time tracers are due to the application of the imaging veto masks, which remain the same as applied to DR1. The area of the BGS sample has grown by 65% compared to DR1, and this is simply due to the footprint of the newly observed tiles.

Tracer	# of good $z$	$z$ range	Area [deg$^{2}$]	$C_{\mathrm{assign}}$	$z$ succ.
BGS	1,188,526	0.1 − 0.4	12,355	75.5%	98.8%
LRG	4,468,483	0.4 − 1.1	10,031	82.6%	99.0%
ELG	6,534,844	0.8 − 1.6	10,352	53.7%	73.9%
QSO	2,062,839	0.8 − 3.5	11,181	93.6%	68.0%

The BGS sample we use is nearly $4\times$ the size of the DR1 BGS sample used for BAO measurements. The increased area and completeness account for half of this increase. The rest of the increase is due to a relaxation of the absolute magnitude cut. First, we changed the definition of the absolute magnitudes used for the cut. In this analysis, we do not apply any $k$ or $e$ corrections and we instead determine the absolute $r$-band magnitudes, $M_r$ purely based on the apparent magnitude of the galaxy and the distance modulus to its redshift.³ Based on this new definition, we also changed the absolute magnitude cut to the value at which no galaxies would be removed from the sample at its maximum redshift of 0.4. This was determined to be $M_r < -21.35$, whereas the cut applied in DR1 was $M_r < -21.5$.

Figure 3 shows the comoving number density of the BGS sample (applying completeness corrections) that we use in solid light green. The sample used in our DR1 analysis is shown in dot-dashed dark green. The total BGS_BRIGHT parent sample is shown with a dashed black curve. For the total sample, the number density quickly becomes greater than the range shown on the plot, because the target sample was chosen with a simple $r$-band flux cut [59]. The density of the sample chosen for this analysis matches that of the parent sample at the upper redshift limit $z = 0.4$. The BGS sample we use has approximately double the number density of the one used in the DR1 analysis. Our lack of use of $k$ or $e$ corrections makes the number density evolve with redshift slightly more than for the DR1 BGS sample. However, the change of number density with redshift is $\sim$25%, which is still lower than that of the other tracers used in the analysis. The DESI DR2 data contains more than 9.1 million BGS galaxies (considering both the BGS_BRIGHT and BGS_ANY samples [59]) with good redshifts $0.1 < z < 0.4$, of which we use only a small fraction (Table II). Using all of the BGS galaxies instead would increase the effective volume by approximately 20%. The exact amount depends on the assumed clustering amplitude, which will change significantly with redshift, as the effective luminosity threshold of the sample decreases at lower redshift. This large change in the physical properties of the sample introduces significant modeling complexities as many elements of our analysis pipeline assume a roughly constant clustering amplitude (e.g., reconstruction, covariance estimation, BAO fitting). We aim to utilize more of the BGS in future DESI cosmological analyses.

Figure Figure 4 shows the raw and the completeness-corrected comoving number density for each of the samples used in this analysis. Dotted vertical lines denote the redshift binning we apply to obtain BAO measurements. The labels, number of redshifts, and effective volumes⁴ for the data in each of these redshift bins are provided in Table III. The redshift binning is the same as applied to DR1 and one can thus compare directly to the values in table 2 of [21]. The effective volume has increased compared to DR1 by a factor that varies from 1.8 (QSO) to 3.1 (ELG2).

Tracer	Redshift range	$N_{\mathrm{tracer}}$	$V_{\mathrm{eff}}$ (Gpc$^3$)
BGS	0.1 − 0.4	1,188,526	3.8
LRG1	0.4 − 0.6	1,052,151	4.9
LRG2	0.6 − 0.8	1,613,562	7.6
LRG3+ELG1	0.8 − 1.1	4,540,343	14.8
ELG2	1.1 − 1.6	3,797,271	8.3
QSO	0.8 − 2.1	1,461,588	2.7
Ly$\alpha$	1.8 − 4.2	1,289,874	—
LRG3	0.8 − 1.1	1,802,770	9.8
ELG1	0.8 − 1.1	2,737,573	5.8

The version of the LSS catalogs used in this analysis, DR2 v1.1/BAO, will be released with DESI DR2. Prior to producing the results presented in this paper, the DESI BAO analysis team further processed the LSS catalogs

to blind the results and apply BAO reconstruction. The blinding was removed only after the analysis choices were validated and fixed. We provide more detail on each process below.

Blinding

To minimize the risk of confirmation bias in our analysis, we used a blinding scheme to deliberately conceal the true position of the BAO peak observed from the data until all choices about our inference pipeline were finalized. This blinding was applied at the level of the LSS catalogs, using the same prescription as was applied to DR1 and detailed in [68], but with a new random seed to produce unknown (but controlled) shifts to galaxy redshifts that conceal the true BAO peak position. A number of validation tests of both the data and the analysis pipeline were performed on the blinded data catalogs, as described in [48]. Once these tests were passed and the pipeline frozen, the blinding was removed and true catalogs were first processed during the DESI winter collaboration meeting in December 2024.

Reconstruction

Non-linear gravitational evolution induces large-scale bulk flows that can degrade the precision and accuracy of BAO measurements. Density field reconstruction [69] attempts to correct for this effect by estimating the displacement field from the observed galaxy distribution, using it to undo the gravitational flow and partially restore the acoustic peak to its linear regime shape. We ‘reconstructed’ the LSS catalogs using the IterativeFFT reconstruction algorithm [70] implemented in pyrecon⁵ and using the fiducial DR1 settings from [71]. This decision was informed by an extensive comparison of different reconstruction algorithms in the context of DESI [72]. The reconstructed catalogs are used for the clustering and associated BAO measurements from galaxies and quasars presented throughout the rest of this paper.

Lyman- forest catalog

The DESI DR2 quasar catalog for the Ly$\alpha$ forest analysis was constructed in a similar manner to the DR1 catalog described in [22]. The catalog contains high-redshift quasars and is built from the outputs of three automatic classifiers that analyze all quasar targets: the Redrock redshift estimator [57] that classifies based on templates, including quasar templates optimized for DESI quasars [73, 74], a Mg II afterburner that searches for broad Mg II emission in quasar candidates classified as galaxies, and the QuasarNet neural network classifier developed by [75] that has been updated [76] for DESI. The catalog also contains any other targets that are classified by Redrock as high-redshift quasars.

We also catalog two types of absorption systems that are masked by our analysis pipeline: Damped Ly$\alpha$ Absorption Systems (DLAs) and Broad Absorption Line (BAL) quasars. DLAs are systems with high ($N_{\mathrm{HI}} > 2 \times10^{20}\ \mathrm{cm}^{-2}$) column densities of neutral hydrogen in the intergalactic medium (IGM). These systems have somewhat higher clustering than typical Ly$\alpha$ forest absorption and have very broad damping wings that can compromise a significant fraction of a forest spectrum. We use three analysis tools to identify DLAs in the DR2 dataset; the performance of these methods are described in [77]. Approximately 20% of the QSOs in the Ly$\alpha$ forest sample are BAL quasars. These absorption systems contaminate the forest by adding absorption that is uncorrelated with the large-scale structure of the IGM. We identify BALs with the method described in [74]. Further details about the quasar catalog for the Ly$\alpha$ forest are described in [63].

The DR2 quasar catalog has over 1.2 million quasars at $z > 1.77$, including over 820,000 quasars at $z > 2.09$. We use Ly$\alpha$ forests from the spectra of quasars at $z > 2.09$, while quasars at $z > 1.77$ are used as discrete tracers in the computation of the cross-correlation with Ly$\alpha$. Both samples are close to a factor of two larger than the DR1 sample, which had 420,000 quasars at $z > 2.09$ and over 700,000 quasars at $z > 1.77$. In addition to the larger sample size, the average signal-to-noise per pixel is higher, because 67% of DR2 quasars have had at least two observations (out of the four planned for the Ly$\alpha$ quasar sample), while only 43% of DR1 quasars had been observed more than once.

BAO Measurements

Clustering from Galaxies and Quasars

Our baseline BAO measurements from galaxies and quasars are derived from the two-point correlation function (2PCF) in redshift space. We measure it using the Landy-Szalay estimator [78] (adapted to post-reconstruction measurements, as in [79]) using pycorr.⁶ which implements a wrapper around a modified version of the Corrfunc pair-counting code [80]. We bin the correlation function in $s$ and $\mu= \cos\theta$, where $s$ is the redshift-space scalar separation and $\theta$ is the angle between the galaxy pair and the line of sight, using $s$-bins of $4\ h^{-1}\mathrm{Mpc}$ width and 200 $\mu$-bins in $[-1,1]$. We then decompose it into Legendre multipoles, focusing on the monopole and quadrupole moments for the analysis.

We estimate the errors of the correlation function measurements using the RascalC semi-analytical covariance matrix code⁷ [81, 82, 83, 84, 85], in the same way as for the DESI DR1 BAO analysis [21]. The computation gives the contributions of the measured (non-linear) two-point correlation function (including the disconnected four-point function) with survey geometry and selection effects, and shot-noise rescaling calibrated with jackknife resampling to account for missing non-Gaussian contributions. This procedure has been validated using DESI DR1 mocks in [86], which also provides a more complete description of the algorithm. In general, covariance estimation for DR2 is less challenging than DR1 (e.g., the footprint is less complex and the samples have higher completeness). The DESI DR2 covariance pipeline is publicly available on GitHub.⁸

The BAO fitting procedure—described in detail in [21]—involves the use of a template of the correlation function multipoles in a fiducial cosmology (which are converted from Fourier space templates of the power spectrum, [87]). During the parameter posterior sampling, the BAO features in these templates are shifted with respect to those seen in the data, and the amount of shifting is regulated by scaling parameters that are varied freely during the fit. The scaling parameters that shift the BAO features along and across the line of sight are related to the cosmology by

where the fid superscript denotes quantities in the fiducial cosmology that is used to convert redshifts to distances and define the power spectrum template. These factors can be re-parameterized into a different basis,

which modulate an isotropic ($\alpha_{\mathrm{iso}}$) or anisotropic ($\alpha_{\mathrm{AP}}$) shifting of the BAO feature. The measurements are less correlated in the latter basis, and we use it as our baseline when sampling the BAO model posterior. We also define the isotropic BAO distance $D_{\mathrm{V}}(z) \equiv(zD_{\mathrm{M}}(z)^2D_{\mathrm{H}}(z))^{1/3}$, for which we report results in later sections.

The templates are also allowed to vary in amplitude and are combined with a set of nuisance parameters that model the broadband shape of the correlation function. The reported BAO constraints are marginalized over these nuisance parameters so that the information coming from the position of the acoustic feature can be isolated. The broadband parametrization involves a piecewise-spline fitting basis that was introduced in the DR1 analysis, described in detail in [87].

The BAO model and inference pipeline are largely the same as used for DR1 and described in [21], with a few small modifications:

• Scale cut ($s_{\min}$): for the DR1 analysis, the minimum scale used in the BAO fits was $s_{\min} = 50\ h^{-1}\mathrm{Mpc}$. Reference [87] showed through tests on mock catalogs that the recovered values of $\alpha_{\mathrm{iso}}$ and $\alpha_{\mathrm{AP}}$ and their errors are very stable

against changes to the minimum scale in the range $50\,h^{-1}\mathrm{Mpc} \lesssim s_{\min} \lesssim80\,h^{-1}\mathrm{Mpc}$. However, when fitting to the blinded DR2 data we found somewhat large $\chi^2$ values when using $s_{\min}=50\,h^{-1}\mathrm{Mpc}$ for some tracers and redshift bins, which improved significantly when changing the scale cuts to $s_{\min}=60\,h^{-1}\mathrm{Mpc}$ (with only very minor effects on the recovered $\alpha$ values). To guard against the possibility that this effect in the blinded data was due to some unknown systematic affecting the 2PCF in the range $50<s<60\,h^{-1}\mathrm{Mpc}$, we chose to set $s_{\min}=60\,h^{-1}\mathrm{Mpc}$ as our baseline for our BAO fits before unblinding.

• 2D fits for ELG1 and QSO: for tracers with a sufficiently large signal-to-noise ratio in the 2PCF, $\alpha_{\mathrm{iso}}$ and $\alpha_{\mathrm{AP}}$ can be simultaneously determined, using anisotropic clustering information by fitting both the monopole and quadrupole moments of the correlation function (a ‘2D fit’). However, when the signal-to-noise of the quadrupole measurement is lower, robust determination of $\alpha_{\mathrm{AP}}$ is harder. In these cases, we only fit the monopole, which depends on $\alpha_{\mathrm{iso}}$, to avoid the complications of including a weak non-Gaussian constraint on $\alpha_{\mathrm{AP}}$. For the DR1 analysis, our baseline analysis for BGS, ELG1, and QSO used only the monopole for BAO fits. In DR2, the signal-to-noise ratio of the ELG1 and QSO samples has increased sufficiently to promote them to 2D fits for both $\alpha_{\mathrm{iso}}$ and $\alpha_{\mathrm{AP}}$, while for BGS we continue to use the monopole only to measure $\alpha_{\mathrm{iso}}$. This decision was made by verifying the stability of the constraints when fitting the blinded data and mock galaxy catalogs [48].⁹

• Systematic errors: the systematic error treatment in the DESI DR2 analysis incorporates refinements to improve the robustness and accuracy of the results. Fiducial cosmology-related systematic errors were adjusted, increasing the systematic uncertainty on $\alpha_{\mathrm{AP}}$ from 0.1% in DR1 to 0.18%, based on including an evolving dark energy model motivated by the DR1 results in the test [88]. Systematics related to unknown details of the small-scale galaxy-halo connection, which dominated the DR1 error budget [89, 90], were refined to be tracer-specific. The contributions for LRG, ELG, QSO, and BGS tracers were updated to better capture effects dependent on redshift and galaxy type. The details of these choices are presented in Section V B of [48]. The systematic error in $\alpha_{\mathrm{iso}}$ increases the total error budget in this parameter by a fractional amount between 1% (for BGS) and 9% (LRG3+ELG1), while for $\alpha_{\mathrm{AP}}$ the increase ranges between 0.1% (QSO) and 2% (LRG3+ELG1). Thus, in all cases, our error budget remains dominated by statistical errors associated with finite survey volume and sampling density. Additionally, we assessed the impact of including theoretically motivated systematic correlations across redshift bins and found no significant effect on our main results.

The updated methods, systematic error estimates and the range of validation tests performed on blinded data before the pipeline was frozen are described in detail in the supporting publication [48].

Figure 5 shows the correlation function multipoles around the BAO scale measured from the reconstructed DESI DR2 catalogs (first eight panels, starting from the upper left). For those data vectors that are used for the cosmological inference presented in the following sections, we also show the best-fit BAO model as a solid line. The acoustic feature is successfully detected as a distinct peak in the correlation around 100 $h^{-1}\mathrm{Mpc}$ for all tracers. The statistical significance of the detection ranges from 5.6$\sigma$ for the QSO, to 14.7$\sigma$ for LRG3+ELG1, which constitutes the strongest detection of the BAO feature from a galaxy survey to date.

Clustering from the Lyman- Forest

The DR2 Ly$\alpha$ forest BAO measurement [63] is based on the combination of four correlation functions: the auto-correlation function of the Ly$\alpha$ forest in region A (1040–1205 Å; Ly$\alpha$(A)$\times$Ly$\alpha$(A)), the auto-correlation of the Ly$\alpha$ forest in region B (920–1020 Å; Ly$\alpha$(A)$\times$Ly$\alpha$(B)), the cross-correlation of the forest in region A with quasars (Ly$\alpha$(A)$\times$QSO), and of region B with quasars (Ly$\alpha$(B)$\times$QSO). We measure the Ly$\alpha$ forest in these two distinct regions because region B has additional astrophysical noise due to higher-order Lyman series absorption, and the spectra of region B are generally lower signal-to-noise ratio because this region is only observed in higher-redshift quasars. We do not measure the auto-correlation of region B nor the auto-correlation of the quasars because these have even lower signal-to-noise ratio. For regions A and B we measure the forest absorption relative to an estimate of the forest continuum with the picca code.¹⁰ This code fits a mean continuum model plus two diversity parameters for each quasar to derive the over or under-density of the forest absorption.

We measure these four correlation functions with the picca analysis code, along with the covariance and cross-covariance of the four correlation functions. The correlations are calculated from 0–200 $h^{-1}\mathrm{Mpc}$ in 4 $h^{-1}\mathrm{Mpc}$ bins in $r_{\parallel}$ and $r_{\perp}$ that correspond to comoving distances along and transverse to the line of sight, respectively. We fit the correlation functions with the Vega package¹¹

measure the two BAO parameters $(\alpha_{\parallel}, \alpha_{\perp})$ and 15 nuisance parameters that include the Ly$\alpha$ forest bias, the Ly$\alpha$ redshift-space distortion (RSD) parameter, quasar bias, the biases of various metals in the IGM, and several others. All four correlation functions are fit in 2D from $30 < r < 180\,h^{-1}\,\mathrm{Mpc}$. The fit has 9306 data points and 17 parameters. Further details are given in [63] and references therein. Figure 5 shows the Ly$\alpha$ autocorrelation function and the cross-correlation with quasars in the bottom right panel. This panel compresses these 2D measurements into a single correlation, while the analysis in [63] uses the whole 2D information without this compression.

We finalized the Ly$\alpha$ forest analysis choices with blinded data and mock datasets before we measured the BAO parameters. Unlike for lower-redshift galaxies and quasars, the blinding procedure for the Ly$\alpha$ forest analysis shifts the BAO peak location in the correlation function (following the method developed for DR1, see [22]). Briefly, we calculated a correlation function model that sets nuisance parameters to the values in DR1 and then computed a blinding template with Vega for some shift $\Delta\alpha_{\parallel}, \Delta\alpha_{\perp}$ that our analysis pipeline (picca) applied to any calculation of the correlation function. The shift values were not stored anywhere, and instead we stored the random number seed and other information necessary to recover the shift. The virtue of blinding the correlation function, rather than the catalog-level blinding used at lower redshifts, is that it preserves the locations of features in the correlation function due to metal-line contamination, which otherwise would have revealed the magnitude and direction of the blinding.

There are two substantial improvements in the Ly$\alpha$ analysis relative to DESI DR1 [22] that we describe in the companion paper [63] and additional, important improvements that we describe in two supporting papers [77, 91]. The two substantial improvements are in how we account for the distortion of the correlation function due to the continuum fitting, and the modeling of metal contamination in the forest. Supporting paper [91] describes improvements to the number and quality of mock catalogs used to validate our analysis pipeline. The most significant improvement in the mocks is the inclusion of non-linear broadening of the BAO peak. Another important improvement is the development of a new software package to identify DLAs based on template fitting, as well as a careful characterization of the purity and completeness of our DLA catalog with mock datasets. We now use the DLA template-fitting package, as well as two other codes based on a Convolutional Neural Network [92] and a Gaussian Process method [93], to produce the DLA catalog. The new template-fitting DLA finder, the performance of the three finders, and the construction of the final DLA catalog are described in another supporting paper [77]. As shown in [63], the cumulative impact of all these improvements on the BAO results is very small.

The Ly$\alpha$ forest measures the BAO signal much better in the line-of-sight direction than in the transverse direction due to the large value of the Ly$\alpha$ RSD parameter. The optimal combination of $\alpha_{\parallel}$ and $\alpha_{\perp}$ for the Ly$\alpha$ dataset is $\alpha_{\parallel}^{0.55}\alpha_{\perp}^{0.45}$, which differs from the best-measured $\alpha_{\mathrm{iso}}$ quantity for the galaxy and quasar BAO in eq. (12); we measure this optimal combination with a precision of 0.64% at an effective redshift of $z_{\mathrm{eff}} = 2.33$ with the DR2 Ly$\alpha$ sample, a significant improvement relative to even the 1.1% measurement from DESI DR1. Nevertheless, in this paper we still report Ly$\alpha$ BAO results in the conventional $\alpha_{\mathrm{iso}}$, $\alpha_{\mathrm{AP}}$ basis as for the other tracers.

Several recent theoretical studies have reported $\sim3\sigma$ measurements of a shift in the BAO position in the Ly$\alpha$ forest due to non-linear evolution [94, 95, 96]. Each of these studies has various strengths and weaknesses, and do not agree on the magnitude and direction of the shift, so we have adopted a theoretical systematic uncertainty of $\Delta\alpha_{\parallel}=0.3\%$, $\Delta\alpha_{\perp}=0.3\%$, which approximately encompasses the theoretical estimates. We add this systematic uncertainty to the statistical covariance matrix to obtain the total error budget. This systematic error reduces the precision of the isotropic BAO measurement to 0.7%. We discuss the theoretical studies and our choice of systematic error further in [63].

DESI Distance Measurements

As implied by eq. (11), once the BAO scaling parameters are fit through a template approach, we can use the assumed fiducial cosmology to convert our measurements into cosmological distance estimates. Table IV presents the distance measurements derived from the DESI DR2 BAO measurements, reported at the effective redshift $z_{\mathrm{eff}}$ (defined in Eq. (2.1) of [21]) for each bin. This definition of $z_{\mathrm{eff}}$ is chosen to best represent the redshift of the maximum statistical weight of each sample. However, as long as the conversion of the measured BAO $\alpha$ parameters to distances at this redshift is correctly performed, small changes to the definition of $z_{\mathrm{eff}}$ are inconsequential for cosmology inference. For BGS, we only report the angle-averaged distance ratio $D_{\mathrm{V}}/r_{\mathrm{d}}$, since we performed a purely isotropic BAO fit to the monopole. For the other tracers, we report the marginalized constraints on angle-averaged distance, $D_{\mathrm{V}}/r_{\mathrm{d}}$, and the anisotropic factor, $D_{\mathrm{M}}/D_{\mathrm{H}}$, that is inversely proportional to $\alpha_{\mathrm{AP}}$. We also convert these into perpendicular and parallel distances. The reported constraints include the contribution from the systematic error budget that is described in § 3, and the final distance measurements and uncertainties do not depend on the choice of fiducial cosmology (see Section VII of [48]).

Tracer	$z_{\mathrm{eff}}$	$\alpha_{\mathrm{iso}}$	$\alpha_{\mathrm{AP}}$	$D_V/r_d$	$D_M/D_H$	$r_{V,M/H}$	$D_M/r_d$	$D_H/r_d$	$r_{M,H}$
BGS	0.295	$0.9857 \pm0.0093$	—	$7.942 \pm0.075$	—	—	—	—	—
LRG1	0.510	$0.9911 \pm0.0077$	$0.9555 \pm0.0261$	$12.720 \pm0.099$	$0.622 \pm0.017$	0.050	$13.588 \pm0.167$	$21.863 \pm0.425$	$-0.459$
LRG2	0.706	$0.9749 \pm0.0067$	$0.9842 \pm0.0227$	$16.050 \pm0.110$	$0.892 \pm0.021$	$-0.018$	$17.351 \pm0.177$	$19.455 \pm0.330$	$-0.404$
LRG3+ELG1	0.934	$0.9886 \pm0.0046$	$1.0237 \pm0.0157$	$19.721 \pm0.091$	$1.223 \pm0.019$	0.056	$21.576 \pm0.152$	$17.641 \pm0.193$	$-0.416$
ELG2	1.321	$0.9911 \pm0.0071$	$1.0257 \pm0.0237$	$24.252 \pm0.174$	$1.948 \pm0.045$	0.202	$27.601 \pm0.318$	$14.176 \pm0.221$	$-0.434$
QSO	1.484	$1.0032 \pm0.0153$	$0.9885 \pm0.0564$	$26.055 \pm0.398$	$2.386 \pm0.136$	0.044	$30.512 \pm0.760$	$12.817 \pm0.516$	$-0.500$
Ly$\alpha$	2.330	$0.9971 \pm0.0082$	$1.0071 \pm0.0216$	$31.267 \pm0.256$	$4.518 \pm0.097$	0.574	$38.988 \pm0.531$	$8.632 \pm0.101$	$-0.431$
LRG3	0.922	$0.9936 \pm0.0053$	$0.9996 \pm0.0172$	$19.656 \pm0.105$	$1.232 \pm0.021$	0.106	$21.648 \pm0.178$	$17.577 \pm0.213$	$-0.406$
ELG1	0.955	$0.9888 \pm0.0091$	$1.0574 \pm0.0291$	$20.008 \pm0.183$	$1.220 \pm0.033$	0.420	$21.707 \pm0.335$	$17.803 \pm0.297$	$-0.462$

We note that any calculations of the degree of consistency or otherwise of DESI BAO results with those of any external dataset mathematically cannot depend on whether they are performed in terms of $(D_{\mathrm{V}}/r_{\mathrm{d}},D_{\mathrm{M}}/D_{\mathrm{H}})$ or $(D_{\mathrm{M}}/r_{\mathrm{d}},D_{\mathrm{H}}/r_{\mathrm{d}})$, or any other rotated basis in this space.

We show a visualization of the cosmological distance estimates in Figure 6, where our measurements have been normalized by the best-fit DESI $\Lambda$CDM predictions. The previous BAO measurements derived from DESI DR1 [21, 22] are included for comparison. This plot highlights the improvement in the parameter precision of the new data release, as well as the inclusion of a new measurement of the BAO distance ratio (and thus of the perpendicular and parallel distances individually) at $z\sim1.5$ from the QSO (green points).

Internal consistency of DESI BAO measurements

As DR1 is a subset of DR2, and the data reduction and analysis pipelines are very similar, we conservatively assume a perfect correlation between the DR1 and DR2 measurements shown in Figure 6 to assess their consistency. We find improvements in the statistical precision between 30% and 50% for DR2. We find that the largest discrepancies between results are $1.56\sigma$ in $\alpha_{\mathrm{iso}}$ (for LRG3+ELG1) and $1.93\sigma$ in $\alpha_{\mathrm{AP}}$ (for LRG2), with the joint uncertainty defined as $\sigma=\sqrt{\sigma_{\mathrm{DR2}}^{2}+\sigma_{\mathrm{DR1}}^{2}-2C\sigma_{\mathrm{DR1}}\sigma_{\mathrm{DR2}}}$, where $C$ is the correlation coefficient between the two samples.¹² The Kolmogorov-Smirnov (KS) statistic for the differences is 0.25 for 12 data points, corresponding to a $p$-value of 0.40 and indicating high consistency between DR1 and DR2.

In order to check the internal consistency of the BAO measurements in different redshift bins, we interpret each within the flat $\Lambda$CDM cosmological model, using the methods described below in § 4.3. In this model, BAO measurements constrain the parameter $\Omega_{\mathrm{m}}$ and the combination $H_{0}r_{\mathrm{d}}$. Figure 7 shows the results for each individual redshift bin and the combination: the posteriors all have significant overlap, indicating no internal inconsistency within a $\Lambda$CDM model. The largest difference in the recovered parameters is $1.8\sigma$ (between LRG1 and LRG3+ELG1). The joint fit to all tracers returns a best-fit $\chi^{2}/\mathrm{dof}=10.2/(13-2)$.

Comparison to SDSS

In [38] we reported a $\sim3\sigma$ difference between the DESI DR1 value of $D_{\mathrm{M}}/r_{\mathrm{d}}$ measured in the LRG2 redshift bin at $z_{\mathrm{eff}}=0.71$ and the corresponding result at $z=0.7$ previously reported by SDSS [41, 97, 98] from the eBOSS LRG sample (although the discrepancy in the two-dimensional $D_{\mathrm{M}}/r_{\mathrm{d}}$–$D_{\mathrm{H}}/r_{\mathrm{d}}$ plane was less significant). The new DR2 result in this redshift bin lies between the two in the $D_{\mathrm{M}}/r_{\mathrm{d}}$–$D_{\mathrm{H}}/r_{\mathrm{d}}$ plane.

To investigate the consistency, we compare measurements assuming that the coefficient of correlation between the power spectrum measurements in the two surveys, defined as

also describes the correlation between BAO measurements. Here $X_i=(n_iP_i)^2/(1+n_iP_i)^2$ for $i=1,2$, and $X_{12}=(n_1n_2P_1P_2)^2(n_1n_2\sqrt{P_1P_2}+n_1n_2P_1P_2)^{-2}$, where $n_1$, $n_2$, and $n_{12}$ respectively are the mean galaxy number densities of DR2 LRG2, the eBOSS LRGs, and the common sample. $P_1$ and $P_2$ are the corresponding power spectrum amplitudes, for simplicity assumed to be the same here. This calculation gives an estimate $C \sim0.57$ for the LRG2 and eBOSS samples—significantly higher than the 0.21 for DR1 due to the larger degree of overlap between the DR2 and SDSS footprints. Assuming this level of correlation between the results, we find the discrepancy between the DR2 and SDSS results has reduced from $3\sigma$ to $\sim2.6\sigma$. Although the assumption of no correlation is less plausible, it sets a lower limit of the discrepancy at the $1.9\sigma$ level (compared to $2.7\sigma$ in DR1). When using the DESI reanalysis of SDSS (as shown in the bottom row of Table 17 in [21]), the discrepancy is $2.3\sigma$ when using $C = 0.57$ and only $1.5\sigma$ when assuming no correlation (compared to $2.8\sigma$ and $2.5\sigma$ for DR1).

For a systematic comparison between DR2 and SDSS, we compare the results in four redshift bins (LRG1, LRG2, QSO and Ly$\alpha$) where the effective redshifts of DESI and SDSS tracers are similar ($\Delta z_{\rm eff} < 0.03$). We convert the SDSS results into $\alpha_{\rm iso}$ and $\alpha_{\rm AP}$ using the reported correlations between $\alpha_\perp$ and $\alpha_\parallel$, and then compute the significance of differences between these values and the DR2 results assuming no correlation between the surveys.

The biggest difference we find is $1.49\sigma$ for $\alpha_{\rm iso}$ and $1.69\sigma$ for $\alpha_{\rm AP}$. The KS statistic is $0.30$, $p = 0.39$. If we conservatively treat SDSS as a pure subsample of DESI DR2 and use the same method to calculate correlation coefficients as in the DR1 vs DR2 comparison above, we find $p = 0.15$ from the KS test. We conclude that there is no significant discrepancy between the DR2 measurements and those from SDSS.

Unlike the case for DR1, in DR2 the effective volume and statistical constraining power of each DESI redshift bin is far larger than that of the corresponding SDSS counterpart. We therefore no longer consider any combination of bins picked from DESI and SDSS data at different redshifts as done in [38].

External Data

We use results from a number of external experiments together with our BAO measurements in order to obtain the most precise cosmological constraints. These are briefly described below.

Big Bang Nucleosynthesis (BBN) prior on $\omega_{\rm b}h^2$

Absent an external determination of $r_{\rm d}$, BAO measurements serve as an uncalibrated ruler and therefore measure the degenerate combination $H_0r_{\rm d}$ rather than $H_0$ and $r_{\rm d}$ individually. Given that BAO also measure $\Omega_{\rm m}$, from eq. (2) we can see that, assuming standard neutrino content ($N_{\rm eff}=3.044$ and $\sum m_\nu=0.06\ {\rm eV}$), knowledge of $\Omega_{\rm b}h^2$ is sufficient to break this degeneracy. Careful observational determinations of the primordial deuterium abundance D/H [99] and the helium fraction $Y_{\rm P}$ [100] in astrophysical systems can be connected to predictions from BBN for the abundances of D and ${}^{4}{\rm He}$ in the early Universe to determine the baryon-to-photon ratio and thus the physical baryon density $\Omega_{\rm b}h^2$, when combined with a measurement of the CMB temperature [32].

As in DR1, we use estimates of $\Omega_{\rm b}h^2$ determined by [101] based on the PRyMordial code [102] including marginalization over uncertainties in nuclear reaction rates. This corresponds to

in the standard $\Lambda$CDM model, and

when allowing $N_{\rm eff}$ to vary (in which case there is also a covariance between $\Omega_{\rm b}h^2$ and $N_{\rm eff}$). We use this prior when we require calibrated results that are independent of CMB information, and denote it as BBN in plots and tables.

Cosmic microwave background (CMB)

The power spectra of anisotropies in the CMB contain a wealth of information on cosmological parameters that is complementary to that from BAO. The large-scale temperature and polarisation power spectra and CMB lensing power spectrum have been exquisitely measured by Planck [24, 103], with additional information from smaller scales and improved lensing reconstruction provided by the Atacama Cosmology Telescope (ACT; [104, 105, 106]) and the South Pole Telescope (SPT; [107]).

When using the full power of the CMB information, our baseline analysis makes use of the temperature ($TT$), polarization ($EE$) and cross ($TE$) power spectra from Planck, specifically using the small, Commander (for $\ell< 30$) and CamSpec (for $\ell\ge30$) likelihoods, plus the combination of Planck and ACT DR6 CMB lensing from [104].¹³ CamSpec [108, 109] is a new likelihood built on the latest NPIPE PR4 data release from the Planck collaboration, which replaces the original Plik likelihood based on the older PR3 release and includes some important differences in methodology. Another set of likelihoods based on PR4, known as LoLLiPoP (for low $\ell$) and HiLLiPoP (for high $\ell$) have also been independently released [110], which for brevity we will refer to as L–H in the following. In [38] [46] we used Plik as our default high-$\ell$ likelihood but noted that both CamSpec and L–H reduced the so-called $A_{\rm L}$ anomaly which has a small effect on certain cosmological parameters. They also make use of slightly more data, with larger sky fractions, than Plik. Our choice of CamSpec as the default is based on the speed of evaluation of the likelihood code relative to L–H, but where any results depend on the choice of likelihood (notably for the sum of neutrino masses $\sum m_\nu$) we provide both sets of results.

Fits using the full CMB likelihoods require specification of a full physical model, including computation of perturbations, the late integrated Sachs–Wolfe (ISW) effect and lensing-induced distortions. However, certain aspects of the early Universe can be robustly constrained independently of any modifications of the late-time cosmological model. The most precisely determined quantity from the CMB is the angular scale of the acoustic fluctuations, $\theta_* = r_*/D_{\rm M}(z_*)$, where $r_*$ is the comoving sound horizon at recombination and $D_{\rm M}(z_*)$ is the transverse comoving distance to that redshift. This quantity is the direct analogue of the measurement of $D_{\rm M}/r_{\rm d}$ that we make using BAO. We can consider ‘BAO-only’ information by combining galaxy BAO with $\theta_*$. When doing this, we use a Gaussian prior on $\theta_*$,

The mean value of this prior matches the Planck result for $\Lambda$CDM, and the width of the prior has been conservatively increased by $\sim75\%$ over the Planck reported uncertainty, to account for the (small) variation seen over different assumed late-time models.

More generally, the CMB determines more information than just $\theta_\ast$ independent of the late-time evolution. Following [111], this can be expressed in terms of correlated Gaussian posteriors on the parameters $(\theta_\ast,\omega_b,\omega_{bc})$ obtained after marginalizing over the ISW and lensing contributions and other possible late-time effects.¹⁴ We use this information, in the form of a correlated Gaussian prior, as an alternative to using the full CMB likelihood that is also more model-independent. In practice we determine the correlation between parameters from a set of early-Universe results provided by [111] based on the CamSpec likelihood¹⁵ and use this to define the Gaussian prior. Numerical details of the implementation are given in Appendix A. We refer to this in the following text and plots as $(\theta_\ast,\omega_b,\omega_{bc})_{\mathrm{CMB}}$ priors.

For constraints on dark energy evolution in particular, these CMB priors contain most of the relevant information from the full CMB likelihoods, in the form of a high-redshift calibration of the low-redshift observations like BAO that directly probe the background evolution. In this form, the CMB information is independent of whether or not dark energy evolves. Such evolution would alter the distance to last scattering, which is already allowed for in the definition of $\theta_\ast$ as the angular size of the horizon; the densities $\omega_b$ and $\omega_c$ are determined via a ratio to the radiation density at $z_\ast$, which is independent of evolution at lower redshifts.

Type Ia supernovae (SNe)

Type Ia SNe are luminous standardizable candles that provide another probe of the expansion history of the Universe, especially useful at low redshifts ($0.01 < z < 0.3$) where BAO measurements are limited by cosmic variance. There are three recent SNe datasets and likelihoods available for cosmology: the Pantheon+ [29, 113], Union3 [30] and Dark Energy Survey Year 5 (DESY5; [31]) samples.

The Pantheon+ and Union3 samples consist of compilations of 1,550 and 2,087 spectroscopically-classified SNe, respectively, drawn from multiple observational programs over the past few decades of observations. While 1,363 of the objects are in common between the two, there are differences in the calibration, modeling and treatment of systematic uncertainties.

The DESY5 sample consists of 1,635 photometrically-classified $z > 0.1$ SNe with uniform calibration drawn from a single survey up to $z = 1.13$, together with a small set of 194 SNe at $z < 0.1$ drawn from historical sources, some of which are in common with the Pantheon+ and Union3 datasets. Despite the larger differences in the underlying data, the Pantheon+ and DESY5 analysis methodologies are quite similar in nature, although DESY5 uses the SALT3 [114] light-curve fitting model as opposed to SALT2 [115], among other differences. Union3 shares a larger fraction of the data with Pantheon+ but uses a different approach based on the Bayesian hierarchical modeling framework Unity1.5 [116].

As in our DR1 analysis, we do not choose between these three options but instead present our headline results in combination with each SNe dataset independently, noting the differences in interpretation where appropriate. All three likelihoods are implemented in the Cobaya sampling code [117] [118] and the underlying data are publicly available.¹⁶

In Section VII, we illustrate the distance modulus residuals for the SNe datasets in redshift bins. For visual purposes only, we rebin the SNe in redshift using the same bins for all datasets, and then calculate the weighted average distance moduli and errors using the inverse of the covariance matrix (including both statistical and systematic errors and after removing the weighted mean of the distance moduli). All cosmological fits use the data and covariance as originally provided. It is important to note that in all SN samples the absolute magnitude $M$ of the SNe is completely degenerate with the Hubble parameter $H_0$ and therefore, unless some external calibration of either $M$ or $H_0$ is assumed, the distance moduli can be scaled by any constant offset. For cosmological inference we therefore always marginalize over the value of $M$; for plots we show the residuals of the distance modulus relative to the mean over all SNe.

Dark Energy Survey 32pt

Weak gravitational lensing of galaxies and its combination with galaxy clustering are local probes of structure formation. For models with an evolving dark energy equation of state, we use galaxy-galaxy, galaxy-shear, and shear-shear two-point correlation functions (referred to here as 3×2pt) measurements from the Dark Energy Survey Year-3 (DESY3) analysis [119]. These data vectors were estimated from observations of over 10 million lens galaxies in the MagLim sample covering an area of approximately 4000 square degrees. We adopt the same modeling choices and scale cuts as the corresponding DESY3 analysis [120]. We choose the same priors as DESY3 with the exception of the parameters $w_0$, $w_a$, $H_0r_d$ and $\Omega_m$, where we impose priors that match those used for the DESI analysis in this paper. In particular, in post-processing we compute $H_0r_d$ as a derived parameter and reweight the chains to impose the same prior for free parameters in the DESI chains. We sample the DESY3 likelihood using the publicly available DESY3 CosmoSIS pipeline [121]. To combine these constraints with DESI and SNe data at the posterior level, we use CombineHarvesterFlow [122] to fit normalizing flows to the DESY3 chains and reweight the DESI and DESI+SNe chains.

Cosmological Inference

We use the cosmological inference code Cobaya [117] [118] to sample posteriors in parameter space for cosmological inference, using Metropolis-Hastings Monte Carlo Markov Chain (MCMC) sampling. Theory models are computed using interfaces to CAMB [123] and Cobaya likelihoods. When fitting to BAO data only, we sample over $\Omega_m$ and the parameter $h r_d$, where $h = H_0/(100\ {\rm km\ s^{-1}\ Mpc^{-1}})$. When the BAO data are calibrated through the use of an external prior such as BBN, we sample in $H_0$ and $\omega_b$ instead of $h r_d$. Runs involving use of the full CMB likelihood sample the base set of six cosmological parameters $(\omega_{bc},\omega_b,100\theta_{\rm MC},\ln(10^{10}A_s),n_s,\tau)$, where $\theta_{\rm MC}$ is an approximation¹⁷ to the acoustic angular scale $\theta_*$ based on fitting formulae given in [124], $\ln(10^{10}A_s)$ and $n_s$ are the amplitude and spectral index of the primordial scalar perturbations, and $\tau$ is the optical depth. Extended cosmological models additionally allow the spatial curvature $\Omega_K$, the dark energy equation of state parameters $w_0$ (or $w$) and $w_a$, or the sum of neutrino masses $\sum m_\nu$ to vary in addition to the base parameters mentioned above. When $\sum m_\nu$ is varied, we assume three degenerate mass eigenstates unless testing specific mass ordering scenarios; when it is not varied, it is fixed to a default value of 0.06 eV assuming a single non-zero mass eigenstate. Prior ranges on all sampled parameters match those given in Table 2 of [38].

The convergence criterion for MCMC sampling is that the Gelman-Rubin statistic [125] satisfies $R-1<0.01$, or that the chains have an effective sample size of $\gtrsim10^3$, whichever is longer. Summary statistics for our chains as well as plots are obtained with getdist¹⁸ software package [126]. For 1D marginalized posterior results we

quote the mean and standard deviation when the distributions are symmetric and the 68% minimal credible interval when they are not. Where only limits on parameter values can be determined, upper or lower bounds are quoted at the 68% level, except for $\sum m_\nu$ where we quote the 95% upper bound to ensure comparability with previous work.

In order to determine the best fit points and the corresponding $\chi^2$ we use the iminuit [127] algorithm starting from the maximum a posteriori (MAP) points of each of the chains in the MCMC sampling. When comparing the fits of two different models we compare the quantity $\Delta\chi^2_{\rm MAP}\equiv-2\Delta\ln\mathcal{L}$ representing twice the difference in the negative log posteriors at the maximum posterior points for each model. This measure accounts for the contribution of any non-uniform priors when evaluating the differences between two MAP points. For likelihood combinations including DESY3 ($3\times2\mathrm{pt}$), when computing the best-fit points and $\chi^2$, we directly use the version of the DESY3 likelihood [119] implemented in the CosmoSIS [121] pipeline. For these combinations with DESY3, when calculating the deviance information criterion values, we separate the calculation of the mean $\langle\chi^2\rangle$ via the sum $\langle\chi^2_{\mathrm{DESI(+SNe)}}\rangle+\langle\chi^2_{\mathrm{DESY3}}\rangle$ and compute each term with the corresponding weighted CombineHarvesterFlow chains. This enables the use of the values of $\chi^2$ from each individual chain, while the distributions being averaged over are equivalent.

Cosmological constraints in the CDM model

We start by presenting cosmological constraints in the base $\Lambda$CDM model, and examining tensions between different datasets in this scenario, before introducing the freedom of extended models in Sections VII and VIII.

As discussed in § 3 C and shown in Figure 7, the values of $\Omega_m$ and $H_0r_d$ inferred from each DESI tracer individually within the $\Lambda$CDM framework are consistent with each other. The result of a combined fit to the BAO results from all redshift bins together is shown in the left panel of Figure 8 compared to those obtained from DR1 and from the CMB. We find

with a correlation coefficient of $r=-0.92$. This represents a $\sim40\%$ improvement in the precision on $\Omega_m$ and $h r_d$ compared to the DR1 results, while being perfectly consistent with them as well as with the SDSS and

DESI+SDSS results reported in [38]. In the following text and figures, we refer to DESI DR2 simply as DESI, unless explicitly comparing to DESI DR1.

Figure 8 shows that the overlap with the CMB posterior has decreased, a sign of the increased discrepancy between the results from DESI BAO and CMB probes when interpreted in base $\Lambda$CDM. We calculate the relative $\chi^{2}$ between the two datasets as

where $\mathbf{p}_{A}$ and $\mathbf{p}_{B}$ are the $(\Omega_{\mathrm{m}}, r_{\mathrm{d}}h)$ values obtained from DESI and the CMB respectively, and $\operatorname{Cov}_{A}$ and $\operatorname{Cov}_{B}$ are the corresponding $2 \times2$ posterior parameter covariances. Converting this $\chi^{2}$ into a probability-to-exceed (PTE) value, we find it is equivalent to a $2.3\sigma$ discrepancy between BAO and CMB in $\Lambda$CDM, increased from $1.9\sigma$ in DR1. However, we note that this reduces to $2.0\sigma$ if CMB lensing is excluded. This discrepancy is part of the reason why more models with a more flexible background expansion history than $\Lambda$CDM, such as the evolving dark energy considered in Section VII below, may provide a statistically preferable fit.

By calibrating the BAO relative distance measurements using external information we are able to determine the Hubble constant $H_0$. The right panel of Figure 8 compares the constraints obtained in the $\Omega_{\mathrm{m}}$-$H_0$ plane from calibrated BAO measurements to those from the CMB. We show DESI BAO results calibrated using the BBN $\omega_{\mathrm{b}}$ prior in the blue contours, while those in orange illustrate the additional combination with the $\theta_{*}$ BAO-like measurement from the CMB. The overlap with the full CMB posterior, shown in pink, is small and has decreased since DR1, showing a $2.3\sigma$ discrepancy with DESI after BBN calibration. The centers of both DESI+BBN and DESI+BBN+$\theta_{*}$ posteriors remain quite well aligned with the degeneracy direction of the CMB, which is $\Omega_{\mathrm{m}}h^{3}$ [128] and shown by the black dotted line in the figure, while the degeneracy direction of the DESI+BBN+$\theta_{*}$ contour instead follows the line of constant $\Omega_{\mathrm{m}}h^{2}$.

Using the conservative BBN prior on $\omega_{\mathrm{b}}$, we obtain

a 0.8% precision measurement that is independent of any CMB anisotropy information and comparable to the result from the CMB [24]. This value is 28% more precise than the corresponding result from DR1, $H_{0} = (68.53 \pm0.80)\ {\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$, but very consistent with it. While the result is clearly more in agreement with the low values obtained from early-Universe measurements than with those from the Cepheid-calibrated local distance ladder [129], it is also noticeably higher than the Planck value, a reflection of the growing tension between DESI and Planck when interpreted in flat $\Lambda$CDM. Adding information on the very well-measured acoustic angular scale further improves this result to

now slightly more precise than from the CMB alone [110]. Figure 9 shows the DESI+BBN result for $\Omega_m$ and $H_0$ relative to the SH0ES result [130]. The contours also show how the constituent tracers of the DR2 sample at different redshifts contribute to the final constraint, with the degeneracy directions of the contours changing as the best measured combination of transverse and line-of-sight BAO changes with redshift. In $\Lambda$CDM the tension between the DESI+BBN and SH0ES $H_0$ results now stands at $4.5\sigma$ independent of the CMB. Note that the DESI+BBN result does assume standard pre-recombination physics to determine $r_d$ through eq. (2).

We have highlighted the tension between DESI and CMB in $\Lambda$CDM in order to provide context to the results for extended models in the following sections. However, given that this tension is not close to $3\sigma$, it is still valid to combine the two datasets within the $\Lambda$CDM model and obtain joint constraints. In this case we find

with a correlation coefficient of $r = -0.975$.

We also allow for spatial curvature to vary in our cosmological fits and we do not find a significance preference for a non-flat $\Lambda$CDM model. Table V summarizes the cosmological parameter results from DESI alone as well as in combination with external datasets, in both $\Lambda$CDM and extended models.

Model/Dataset	$\Omega_{\mathrm{m}}$	$H_0$ [km s$^{-1}$ Mpc$^{-1}$]	$10^3\Omega_{\mathrm{K}}$	$w$ or $w_0$	$w_a$
$\Lambda$CDM
CMB	$0.3169 \pm0.0065$	$67.14 \pm0.47$	—	—	—
DESI	$0.2975 \pm0.0086$	—	—	—	—
DESI+BBN	$0.2977 \pm0.0086$	$68.51 \pm0.58$	—	—	—
DESI+BBN+$\theta_*$	$0.2967 \pm0.0045$	$68.45 \pm0.47$	—	—	—
DESI+CMB	$0.3027 \pm0.0036$	$68.17 \pm0.28$	—	—	—
$\Lambda$CDM+$\Omega_{\mathrm{K}}$
CMB	$0.354^{+0.020}_{-0.023}$	$63.3 \pm2.1$	$-10.7^{+6.4}_{-5.3}$	—	—
DESI	$0.293 \pm0.012$	—	$25 \pm41$	—	—
DESI+CMB	$0.3034 \pm0.0037$	$68.50 \pm0.33$	$2.3 \pm1.1$	—	—
$w$CDM
CMB	$0.203^{+0.017}_{-0.060}$	$85^{+10}_{-6}$	—	$-1.55^{+0.17}_{-0.37}$	—
DESI	$0.2969 \pm0.0089$	—	—	$-0.916 \pm0.078$	—
DESI+Pantheon+	$0.2976 \pm0.0087$	—	—	$-0.914 \pm0.040$	—
DESI+Union3	$0.2973 \pm0.0091$	—	—	$-0.866 \pm0.052$	—
DESI+DESY5	$0.2977 \pm0.0091$	—	—	$-0.872 \pm0.039$	—
DESI+CMB	$0.2927 \pm0.0073$	$69.51 \pm0.92$	—	$-1.055 \pm0.036$	—
DESI+CMB+Pantheon+	$0.3047 \pm0.0051$	$67.97 \pm0.57$	—	$-0.995 \pm0.023$	—
DESI+CMB+Union3	$0.3044 \pm0.0059$	$68.01 \pm0.68$	—	$-0.997 \pm0.027$	—
DESI+CMB+DESY5	$0.3098 \pm0.0050$	$67.34 \pm0.54$	—	$-0.971 \pm0.021$	—
$w_0w_a$CDM
CMB	$0.220^{+0.019}_{-0.078}$	$83^{+20}_{-6}$	—	$-1.23^{+0.44}_{-0.61}$	$< -0.504$
DESI	$0.352^{+0.041}_{-0.018}$	—	—	$-0.48^{+0.35}_{-0.17}$	$< -1.34$
DESI+Pantheon+	$0.298^{+0.025}_{-0.011}$	—	—	$-0.888^{+0.055}_{-0.064}$	$-0.17 \pm0.46$
DESI+Union3	$0.328^{+0.019}_{-0.014}$	—	—	$-0.70 \pm0.11$	$-0.99 \pm0.57$
DESI+DESY5	$0.319^{+0.017}_{-0.011}$	—	—	$-0.781^{+0.067}_{-0.076}$	$-0.72 \pm0.47$
DESI+$\left(\theta_*,\omega_b,\omega_{bc}\right)_{\mathrm{CMB}}$	$0.353 \pm0.022$	$63.7^{+1.7}_{-2.2}$	—	$-0.43 \pm0.22$	$-1.72 \pm0.64$
DESI+CMB (no lensing)	$0.352 \pm0.021$	$63.7^{+1.7}_{-2.1}$	—	$-0.43 \pm0.21$	$-1.70 \pm0.60$
DESI+CMB	$0.353 \pm0.021$	$63.6^{+1.6}_{-2.1}$	—	$-0.42 \pm0.21$	$-1.75 \pm0.58$
DESI+CMB+Pantheon+	$0.3114 \pm0.0057$	$67.51 \pm0.59$	—	$-0.838 \pm0.055$	$-0.62^{+0.22}_{-0.19}$
DESI+CMB+Union3	$0.3275 \pm0.0086$	$65.91 \pm0.84$	—	$-0.667 \pm0.088$	$-1.09^{+0.31}_{-0.27}$
DESI+CMB+DESY5	$0.3191 \pm0.0056$	$66.74 \pm0.56$	—	$-0.752 \pm0.057$	$-0.86^{+0.23}_{-0.20}$
DESI+DESY3 (3$\times$2pt)+Pantheon+	$0.3140 \pm0.0091$	—	—	$-0.870 \pm0.061$	$-0.46^{+0.33}_{-0.29}$
DESI+DESY3 (3$\times$2pt)+Union3	$0.333 \pm0.012$	—	—	$-0.68 \pm0.11$	$-1.09^{+0.48}_{-0.39}$
DESI+DESY3 (3$\times$2pt)+DESY5	$0.3239 \pm0.0092$	—	—	$-0.771 \pm0.068$	$-0.82^{+0.38}_{-0.32}$
$w_0w_a$CDM+$\Omega_{\mathrm{K}}$
DESI	$0.357^{+0.041}_{-0.030}$	—	$-2 \pm56$	$-0.45^{+0.33}_{-0.17}$	$< -1.43$
DESI+CMB+Pantheon+	$0.3117 \pm0.0056$	$67.62 \pm0.60$	$1.1 \pm1.3$	$-0.853 \pm0.057$	$-0.54 \pm0.22$
DESI+CMB+Union3	$0.3273 \pm0.0086$	$65.98 \pm0.86$	$0.6 \pm1.3$	$-0.678 \pm0.092$	$-1.03^{+0.33}_{-0.29}$
DESI+CMB+DESY5	$0.3193 \pm0.0056$	$66.82 \pm0.58$	$0.8 \pm1.3$	$-0.762 \pm0.060$	$-0.81 \pm0.24$

Finally, as in [38], we note a mild to moderate discrepancy between the recovered values of $\Omega_m$ from DESI and SNe in the context of the $\Lambda$CDM model. This is shown in the marginalized posteriors in Figure 10: the discrepancy is $1.7\sigma$ for Pantheon+, $2.1\sigma$ for Union3, and $2.9\sigma$ for DESY5, with all SNe samples preferring higher values of $\Omega_m$ though with larger uncertainties. For $\Lambda$CDM we do not report joint constraints on parameters from any combination of DESI and SNe data. However, as with the CMB, these apparent parameter differences potentially indicate that DESI and at least some of the SNe datasets cannot be consistently fit except with models that have greater freedom in the background evolution, as described in the next section (see also [131]).

Dark Energy

Probing the behavior and nature of dark energy is the primary goal of DESI. The question of perhaps greatest interest, and the one that BAO measurements can best illuminate, is the value of the equation-of-state parameter $w = P/(\rho c^2)$, and its possible evolution with time. To examine this we will primarily use the so-called Chevallier-Polarski-Linder (CPL) parametrization [35, 36] of eq. (9). While this form of $w(a)$ does not arise directly from an underlying physical model, it is a flexible parametrization that is capable of matching the predictions for observable quantities obtained in a wide range of models that are physically motivated [132]. The accompanying paper [49] explores various other parametrizations of $w(z)$, as well as non-parametric reconstruction methods.

For certain ranges of parameters $w_0$ and $w_a$, the parametrization of eq. (9) allows so-called ‘phantom’ behavior of dark energy, in which the equation of state crosses to the regime $w(z) < -1$ [133] where the null energy condition (NEC)—which requires that the energy density of dark energy not increase with the expansion of the Universe—is violated. For single scalar-field models of dark energy, this phantom crossing presents severe theoretical difficulties [e.g., 134, 135]. However, more complex models of dark energy, with multiple fields, other dark energy internal degrees of freedom, or non-minimal coupling, can evade these difficulties, as can some modified gravity models, see, e.g., [136, 137, 138, 139, 140]. We therefore adopt wide uniform priors on the parameters, $w_0 \in\mathcal{U}[-3,1]$ and $w_a \in\mathcal{U}[-3,2]$, together with imposing the condition $w_0+w_a < 0$ to enforce early matter domination. While other justifiable choices are possible, and the values of Bayesian quantities such as the model evidence will always depend on the particular choice used, we consider this the minimal empirical approach. Whenever the equation of state crosses the $w=-1$ boundary we use the parametrized post-Friedmann (PPF) approach of [141, 142] to include dark energy perturbations when calculating CMB power spectra—however, as shown below, the method of accounting for dark energy perturbations does not play a major role, since simply applying an early-Universe CMB prior on $(\theta_*,\omega_b,\omega_c)$ largely reproduces the same results on $w_0$ and $w_a$.

Our primary measure of the statistical significance of preference for evolving dark energy from a given data combination is based on $\Delta\chi^2_{\mathrm{MAP}}$ between the best-fit $\Lambda$CDM and $w_0w_a$CDM models for that combination. Because $\Lambda$CDM is nested within $w_0w_a$CDM, corresponding to $w_0=-1$, $w_a=0$, Wilks’ theorem [143] implies that $\Delta\chi^2_{\mathrm{MAP}}$ should follow a $\chi^2$ distribution with two degrees of freedom under the assumption the null hypothesis ($\Lambda$CDM model) holds, and assuming that errors are Gaussian and correctly estimated. To translate $\Delta\chi^2_{\mathrm{MAP}}$ into familiar terms, we quote the corresponding frequentist significance $N\sigma$ for a 1D Gaussian distribution,

where the left hand side denotes the cumulative distribution of $\chi^2$. We also compute the Deviance Information Criterion (DIC) [144, 145, 146, 147], which takes into account the Bayesian complexity of the model and penalizes including extra parameters.

Results

From DESI DR2 BAO alone, we obtain rather weak constraints on the parameters

which mildly favor the $w_0>-1$, $w_a<0$ quadrant. The upper bound on $w_a$ here is the 68% limit, and $w_a=0$

is not excluded at 95%. As was the case in DR1, BAO data alone define a degeneracy direction in the $w_0$–$w_a$ plane, but they do not show a strong preference for dark energy evolution: the improvement in $\chi^2_{\mathrm{MAP}}$ relative to the $\Lambda$CDM case of $w_0=-1$, $w_a=0$ is equivalent to a preference of just $1.7\sigma$. Note that the posteriors in this poorly constrained case are cut off by the priors, so the marginalized means and limits quoted above are prior-dependent.

The minimal extension we consider, beyond BAO data alone, is to add a high-redshift constraint from the early universe. This can be achieved by imposing CMB-derived priors on $\theta_*$, $\omega_b$ and $\omega_{bc}$, as described in Section IV. These priors are independent of the late-time dark energy, and also marginalize over contributions such as the late ISW effect and CMB lensing. Therefore, they provide us with an early time physics prior that can help us set the sound horizon and is based solely on early-Universe information. The result from this data combination is

While this is still bounded by the $w_a>-3$ prior at the lower end, the posterior already clearly disfavors $\Lambda$CDM. The $\Delta\chi^2_{\mathrm{MAP}}$ value decreases to $-8.0$, indicating a preference for an evolving dark energy equation of state at the $2.4\sigma$ level. Replacing these minimal early-Universe priors with the full CMB information leads to only a small shift in the marginalized posteriors

showing that most of the information that the CMB provides on $w(z)$ comes from its role in anchoring early-Universe values of $(\theta_*,\omega_b,\omega_{bc})$ and thus limiting the freedom for models to fit the low-redshift data without an evolving dark energy component. Nevertheless, when including the full CMB information the $\Delta\chi^2_{\mathrm{MAP}}$ decreases to $-12.5$, corresponding to a $3.1\sigma$ preference for evolving dark energy. This change in the $\Delta\chi^2_{\mathrm{MAP}}$ is driven primarily by the inclusion of CMB lensing, the effect of which is (by construction) not captured in the minimal early-Universe priors (see Appendix A for further discussion and a comparison of posteriors with different choices of CMB likelihoods).

SNe data alone provide a complementary degeneracy direction in the $w_0$-$w_a$ plane, as they measure $w_0$ well independently of $w_a$, which is only weakly constrained. The combination of SNe data with DESI BAO can therefore measure $w_0$ and $w_a$ without having the posteriors cut off by the prior ranges we assumed. The marginalized posterior results are listed in Table V and depend on the choice of SNe dataset, with the significances of the preference for the model over $\Lambda$CDM ranging from $1.7\sigma$ to $3.3\sigma$ as summarized in Table VI.

However, as discussed in [38], the posterior for the combination of DESI BAO and SNe alone allows quite a wide range of posterior values of $\omega_m$ (or, equivalently, $\omega_c$). CMB information places extremely tight constraints on $\omega_m$ that are largely independent of the late-time background model. Therefore, the full statistical power of the data is achieved through the combination of the BAO, CMB and SNe datasets, giving the marginalized posterior results

for the combination with the Pantheon+,

with Union3, and

Datasets	$\Delta\chi^2_{\mathrm{MAP}}$	Significance	$\Delta(\mathrm{DIC})$
DESI	$-4.7$	$1.7\sigma$	$-0.8$
DESI+$(\theta_*,\omega_b,\omega_{bc})_{\mathrm{CMB}}$	$-8.0$	$2.4\sigma$	$-4.4$
DESI+CMB (no lensing)	$-9.7$	$2.7\sigma$	$-5.9$
DESI+CMB	$-12.5$	$3.1\sigma$	$-8.7$
DESI+Pantheon+	$-4.9$	$1.7\sigma$	$-0.7$
DESI+Union3	$-10.1$	$2.7\sigma$	$-6.0$
DESI+DESY5	$-13.6$	$3.3\sigma$	$-9.3$
DESI+DESY3 (3$\times$2pt)	$-7.3$	$2.2\sigma$	$-2.8$
DESI+DESY3 (3$\times$2pt)+DESY5	$-13.8$	$3.3\sigma$	$-9.1$
DESI+CMB+Pantheon+	$-10.7$	$2.8\sigma$	$-6.8$
DESI+CMB+Union3	$-17.4$	$3.8\sigma$	$-13.5$
DESI+CMB+DESY5	$-21.0$	$4.2\sigma$	$-17.2$

with DESY5. The posteriors in these three cases, along with the DESI+CMB posterior, are shown in Figure 11. The $\Delta\chi^2_{\mathrm{MAP}}$ values are $-10.7$, $-17.4$, and $-21.0$, corresponding to preferences for the $w_0w_a\mathrm{CDM}$ model over $\Lambda$CDM at the $2.8\sigma$, $3.8\sigma$, and $4.2\sigma$ levels, for combination with Pantheon+, Union3 and DESY5 respectively. These significances have all increased compared to the values reported in [38] based on the DESI DR1 BAO results.

The deviance information criterion (DIC) values for the combination of DESI+CMB with Pantheon+, Union3 and DESY5 SNe are $\Delta(\mathrm{DIC}) \equiv\mathrm{DIC}_{w_0w_a\mathrm{CDM}}-\mathrm{DIC}_{\Lambda\mathrm{CDM}}=-6.8$, $-13.5$, $-17.2$, respectively. These indicate preferences for the $w_0w_a\mathrm{CDM}$ model consistent with those obtained from the $\Delta\chi^2_{\mathrm{MAP}}$ values above. Again, the changes in the $\Delta(\mathrm{DIC})$ values obtained here with DESI DR2 BAO data compared to the DR1 values reported in [38] show that the preference for $w_0w_a\mathrm{CDM}$ has increased with the additional data. Further details on the calculation of DIC values are given in [49].

The pivot redshift $z_p$ at which $w(z)$ in this parametrization is best constrained by the data depends on the particular combination of datasets used. For DESI+CMB, $z_p=0.53$ and $w_p=w(z_p)=-1.024\pm0.043$: this is a lower pivot redshift and a tighter constraint on $w_p$ than that found for the same combination with DR1 BAO in [38], reflecting the additional constraining power of the DR2 BAO results. For the DESI+CMB+DESY5 combination, we find $z_p=0.31$ and $w_p=-0.954\pm0.024$, indicating a mild preference for a deviation from $w=-1$ at the best-measured redshift. As the other two SNe datasets are slightly less constraining than DESY5, the pivot redshifts for combinations with them are slightly larger, as are the uncertainties on $w_p$, but the results for all choices of DESI+CMB+SNe are mutually consistent.

These results sharpen the preference already seen in [38] for an evolving equation of state for dark energy: although the statistical significances from all data depend somewhat on the choice of SNe dataset included, even the weakest of them (the DESI+CMB+Pantheon+ combination) is still nearly $3\sigma$, and the significance is $3.1\sigma$ even when excluding all SNe data altogether (DESI+CMB). In all cases, the favored $w(z)$ shows a phase of $w>-1$ at low redshifts and a phantom crossing to $w<-1$ above redshifts $z\simeq0.4$. Within the $w_0w_a$CDM model, the necessity of such a crossing and the redshift at which it occurs is determined by the requirement to match the precise CMB measurement of $\theta_*$ together with $\Omega_m$ [148]. The details of the recovered form of $w(z)$ and the $w_0$, $w_a$ parameter values naturally depend on the choice of parametrization. The accompanying paper [49] explores various other parametrizations of $w(z)$ beyond eq. (9), and non-parametric reconstruction methods, that exhibit a similar behavior. Reference [49] also performs binned reconstruction of $w(z)$ without assuming a functional form for the equation of state and finds a consistent picture, as shown in Figure 12. The lowest redshift bin shown in the figure favors a value of $w>-1$ at high significance [49], inconsistent with the $\Lambda$CDM expectation $w=-1$.

The apparent preference for phantom crossing to $w(z)<-1$ at intermediate redshifts, and the consequent violation of the NEC, is thus rather independent of parametrization choices made in the analysis (see, e.g., [149] using DR1 data for non-parametric results). However, the equation of state is not directly observable, and we only observe quantities, such as distances, that depend indirectly on $w(z)$. In some circumstances it may therefore be possible to construct particular models that provide reasonable fits to the low-redshift data while still respecting $w(z)\geq-1$ at all $z$ ([150, 151] provide explicit examples, but [152] discusses the general limitations of such models in fitting DESI+CMB data). The supporting paper [49] examines several NEC-respecting models of dark energy evolution and finds that while they can somewhat outperform $\Lambda$CDM, they have low $\Delta(\mathrm{DIC})$ values relative to $w_0w_a$CDM for the combination of DESI, CMB and SNe data, indicating a preference for phantom crossing.

Table V gives a more complete list of the results for other parameters when fitting this model to different combinations of data. It is worth noting that allowing $w(z)$ to vary does not help resolve the so-called Hubble tension [153], since in $w_0w_a$CDM the recovered $H_0$ is lower than the Planck $\Lambda$CDM value. Although not discussed here, [46] showed that allowing evolving dark energy also does not affect the value of the $S_8$ parameter determined from DESI data, which remains consistent with values from the CMB.

The nature of the evidence for evolving dark energy

The Hubble diagrams in Figure 13 illustrate the nature of the evidence for evolving dark energy and its dependence on the adopted datasets. The upper panels show the isotropic, perpendicular, and parallel BAO measurements ($\alpha_{\mathrm{iso}}$, $\alpha_\perp$, and $\alpha_\parallel$), which are normalized to the predictions of the Planck $\Lambda$CDM cosmology. The lower panels plot $\mu-\mu_{\mathrm{fid}}$, the distance modulus relative to the fiducial Planck $\Lambda$CDM prediction, for the three SNe datasets (noting again that for uncalibrated SNe, the distance modulus is only known to an arbitrary constant offset). Our procedure for creating binned data points from the SNe data is described in Section IV C. Because the fiducial SNe absolute magnitude is unknown, all data points are free to move up or down together by the same amount in $\mu$, and we have chosen the normalization such that error-weighted mean of $\mu-\mu_{\mathrm{fid}}$ is equal to zero. Equivalently, any model curve in these panels can be shifted up or down by a constant $\Delta\mu$, and we have normalized them to match the weighted mean of the data.

In all panels, the horizontal black line represents the prediction of the $\Lambda$CDM model best fit by the CMB, with the $1\sigma$ range of these predictions shown by the shaded gray region. This differs very slightly from the fiducial cosmology based on the Planck 2018 results, due to the switch from Plik to CamSpec and the addition of ACT lensing data, explaining the slight offset from unity. The blue solid curve shows predictions for the $\Lambda$CDM model that best fits the DESI BAO data, while purple represents the model fit to DESI+CMB, which closely matches the best fit to DESI with just the $(\theta_*,\omega_b,\omega_{bc})$ early-Universe priors from the CMB (see Table V). The orange curves in each of the lower panels represent the $\Lambda$CDM models that, from left to right, best fit the DESY5, Union3 and Pantheon+ SNe data respectively. The orange curves in the top row are all those for the $\Lambda$CDM DESY5 best fit.

While the statistical significance of disagreement cannot be judged accurately from this plot alone, the DESI measurements clearly prefer lower distances (by 1-2%) than the Planck $\Lambda$CDM prediction at redshifts $z \le1$. There is a $\Lambda$CDM model that fits the DESI data well (blue curve), but it has a lower $\Omega_m$ than the Planck model (0.297 vs. 0.317), as shown previously in Figure 8. The joint-fit model (purple curve) has an intermediate $\Omega_m = 0.303$, and consequently has a worse fit to both DESI BAO and the CMB (top panels) and also fails to describe the SNe data (lower panels). Similarly, the DESY5 data in the lower left panel exhibit a tension with Planck $\Lambda$CDM, primarily because of the contrast between the low redshift ($z < 0.1$) data at $\mu-\mu_{\rm fid} \approx0.03$ and the points at higher redshift. The story is similar for Union3, but for Pantheon+ the value of $\mu-\mu_{\rm fid}$ at low redshift is smaller, only 0.01.

As is the case for BAO, there exist $\Lambda$CDM models that can reasonably fit the SNe data (orange curves), but these have large $\Omega_m$ that do not well match CMB constraints and are also strongly inconsistent with DESI data, which prefer an $\Omega_m$ value that is lower than Planck, not higher. Conversely, the conflict with the DESIconstrained $\Lambda$CDM models (purple and blue curves) and DESY5 SNe is 0.04 in $\mu$ at low redshifts, $z<0.1$, which is worse than the 0.03 offset from the Planck fiducial model. All these observations point to the fact that the $\Lambda$CDM model struggles to consistently fit all three datasets: BAO, CMB, and SNe.

One way to address this tension is to adopt a model with more flexibility in the background expansion. However, the $w$CDM model with a constant equation of state lacks sufficient flexibility. The cyan dotted curves in the top row and bottom left show a model with constant $w=-0.971$ and $\Omega_m=0.310$, which is the $w$CDM model that best fits the DESI+CMB+DESY5 data combination; in the remaining two panels the cyan curves show the predictions of very similar models obtained by substituting Union3 or Pantheon+ for DESY5. The high-$z$ anchor does not allow enough redshift evolution for this model to provide a good fit to the BAO and SNe data at low redshifts, performing only slightly better than Planck $\Lambda$CDM. If instead the $w$CDM model were chosen to fit the DESI+DESY5 data, it would necessarily have an $\Omega_m$ value that would fail to match the CMB constraints on $(\theta_\ast,\omega_b,\omega_{bc})$.

On the other hand, the $w_0w_a$CDM model does have sufficient flexibility to simultaneously achieve good fits to all three datasets. The green dashed and red dashed curves in the top panels show predictions for $w_0w_a$CDM models with parameters matching the best fits to DESI+DESY5 and DESI+CMB+DESY5 respectively. These are barely distinguishable in the plots, showing that the $w_0w_a$CDM model that best fits the

BAO and SNe automatically also provides a good fit to the CMB. Over the range of redshifts covered by DESI BAO, these model predictions provide a better fit to BAO data than the best such $\Lambda$CDM model—in particular, in the fit to the parallel BAO distances $\alpha_\parallel$—while also simultaneously resolving the mismatch in $\Omega_m$ between DESI and CMB in the $\Lambda$CDM framework. The bottom left panel also shows that they are also equally good at fitting the $\mu$ offsets between low-$z$ and high-$z$ DESY5 data. These observations qualitatively help to explain why these models are strongly preferred over $\Lambda$CDM.

In the equivalent curves in the center and right panels on the bottom row, DESY5 is replaced in all fits by Union3 and Pantheon+, respectively. While the qualitative picture is the same for Union3 as DESY5, the value of $\mu-\mu_{\rm fid}$ at low redshift is smaller for Pantheon+ and so using this dataset does not strengthen the preference for evolving dark energy relative to that already provided by the joint fit to DESI and the CMB.

Are alternative explanations possible?

Given the surprising results from our analysis—our best-fit evolving dark energy models apparently imply phantom dark energy at some redshifts, and order unity changes in $w$ over the redshift range spanned by the observations considered—it is worth considering alternative explanations for the data.

The simplest of these alternatives is simply a statistical fluctuation. Based on the $\Delta\chi^2_{\mathrm{MAP}}$ values presented in Section VII A above, the statistical preference for the best-fit $w_0w_a\mathrm{CDM}$ model over $\Lambda\mathrm{CDM}$ ranges from around $3\sigma$ to over $4\sigma$, depending on the combination of datasets considered. This level of significance is not trivially dismissed, even if it does not yet rise to the $5\sigma$ threshold commonly accepted for establishing new physics. The variations in the statistical significance within this range with the choice of SNe sample also highlight the importance of calibrating the low-to-intermediate redshift SNe distance scale.

The constraining power of SNe in measuring the equation of state comes primarily from the comparison of low-redshift ($z < 0.1$) and high redshift ($z > 0.1$) supernovae. For supernovae at $z > 0.1$, which partially overlap the redshift range of DESI, the $\Lambda\mathrm{CDM}$ model that best fits the DESI data is also a good fit to the SNe data. Relative to models that best fit each of the DESY5, Union3 and Pantheon+ SNe samples alone, over the full redshift range, the DESI best-fit model gives only small shifts in the quality of the fit to the SNe data, with $\Delta\chi^2 = -1.2$, $1.5$ and $2.3$ respectively. Unfortunately, no SNe compilation yet exists in which objects from both redshift regimes were collected as part of a uniform observational program. The Pantheon+ and Union3 datasets are compilations of SNe observed in many individual datasets, each with different calibrations and selection functions. Even for the DESY5 dataset, while 1635 SNe with $z > 0.1$ come from the DES survey program with homogeneous calibration, 194 SNe at $z < 0.1$ are drawn from a mix of historical observational programs with the best-controlled calibrations. The consistent calibration and processing of these data is therefore a delicate operation. Recently, [154] noted an offset in the differences between the standardized brightnesses of low-$z$ and high-$z$ SNe for objects in common between the DESY5 and Pantheon+ samples; however, [155] explain the causes for this and argue that the differences are well justified.

If we exclude the low-$z$ sample entirely and take only the DESY5 SNe, which is the most uniformly calibrated sample, naturally the constraining power and thus the statistical significance of the preference for evolution is reduced, but the best-fit values of $w_0$ and $w_a$ remain far from $\Lambda\mathrm{CDM}$ (central panel of Figure 14). Even if SNe at all redshifts are excluded altogether, the statistical significance of the preference from DESI+CMB alone still exceeds $3\sigma$. Future cosmology analyses of homogeneous SNe samples from ZTF [156, 157] and LSST [158, 159] may shed further light on the relative calibration across redshifts.

Using the simple early-Universe prior on $(\theta_*,\omega_b,\omega_{bc})$ that marginalizes over information dependent on late-time models in addition to DESI gives very similar posterior constraints on $w_0$ and $w_a$ to the full DESI+CMB combination (left panel of Figure 14), and the central $w_0$, $w_a$ values are very stable. However, the tension with $\Lambda\mathrm{CDM}$ drops somewhat (from $3.1\sigma$ to $2.4\sigma$) when using these simple priors instead of the full CMB: this is because the $(\mathrm{Planck}+\mathrm{ACT})$ CMB lensing likelihood makes a sizable contribution to the $\Delta\chi^2_{\mathrm{MAP}}$ calculation.

Although recent independent results from SPT [160] and ACT [104] find cosmological parameter values that are consistent with those from Planck, another possibility is that the CMB constraints used here suffer from some unknown systematic error. It is therefore interesting to consider the constraints that can be obtained entirely independent of the high-redshift CMB anchor. As argued in Section VII B, the primary role of this CMB anchor is that it limits the freedom to vary $\Omega_m$ when fitting to the low-redshift BAO and SNe distance-redshift measurements. This can be achieved instead by replacing the CMB with the DESY3 ($3\times2\mathrm{pt}$) information to obtain a constraint coming entirely from low-redshift cosmological probes. As shown in the right panel of Figure 14, this combination favors the same region of parameter space, although with somewhat larger uncertainties.

We also consider the possibility of an undetected systematic error in our BAO measurements that would have evaded the many internal checks we have performed (cf. details in [21, 48, 63]). For instance, a coherent systematic shift of $1.5\%$ to all the measured $\alpha_{\mathrm{iso}}$ values (eq. (11)), applied in a redshift-independent manner to all tracers, would decrease the tension between DESI and Planck in the $\Lambda\mathrm{CDM}$ model, and thus decrease the evidence for evolving dark energy.²⁰ However, the required magnitude of such a hypothetical shift exceeds the total systematic error budget indicated by our tests by a factor of 6 (for QSO) to 10 (for BGS and LRGs), and we consider this highly unlikely.

Physical models that can satisfactorily explain the data without requiring dark energy evolution are more challenging. Table V confirms the finding, previously reported in [38], that allowing the curvature $\Omega_K$ to vary freely, while increasing uncertainties, does not shift the posterior in $w_0$ and $w_a$ towards $\Lambda\mathrm{CDM}$. While we have fixed $\sum m_\nu=0.06$ eV for our fiducial analysis, it is allowed to vary in Section VIII below, and we find that this has almost no effect on the dark energy constraints. This is because fitting the CMB power spectrum imposes a positive correlation between $\sum m_\nu$ and $\Omega_m$, thus given a physical prior that $\sum m_\nu>0$ eV, in most of the available parameter space neutrino masses can only increase the CMB value of $\Omega_m$, which exacerbates rather than relieves the tension with DESI BAO, and in the region where $0 \leq\sum m_\nu\leq0.06$ eV they have very limited effect. However, the more exotic scenario where $\sum m_\nu$ is treated as an “effective” neutrino mass such that values $\sum m_{\nu,\mathrm{eff}} < 0$ eV are allowed [161, 162, 163] can allow CMB fits to be consistent with lower $\Omega_m$. This has the effect of improving the $\Lambda$CDM fit and thus reducing the preference for evolving dark energy from DESI+CMB+SNe, but not below the level that is found from the relevant DESI+SNe combination alone without any CMB information (which can be up to $3.3\sigma$, Table VI), as shown in Section VIII. That is, allowing negative effective neutrino masses still does not fully remove the preference for non-cosmological constant dark energy. At the same time, such dark energy models do permit $\sum m_{\nu,\mathrm{eff}}>0$ values that are consistent with terrestrial experiments. More standard extensions, such as varying the number of effective neutrino species, $N_{\mathrm{eff}}$, have no effect on the dark energy results.

Arguably, the least contorted way to reconcile our findings with $\Lambda$CDM is to exclude all constraints from SNe data (on the grounds that the significance depends on the choice of SNe dataset and we have no a priori reason to choose one over the others), and then dismiss the remaining $3.1\sigma$ discrepancy with DESI+CMB (or $2.7\sigma$ when also not including CMB lensing) as a statistical fluke.

Neutrinos

Cosmological observations are sensitive to the sum of the neutrino masses. However, this measurement depends on the choice of a cosmological model. So far, our baseline analysis has assumed the sum of neutrino masses to be fixed to the minimum value allowed by terrestrial experiments under a normal mass ordering, $\sum m_\nu=0.06$ eV, modeled as a single massive and two massless states. In this section we allow $\sum m_\nu$ to vary. When doing so, we assume by default three degenerate mass eigenstates except when testing specific physical mass ordering scenarios. This model produces a good approximation of the observable effects of both normal and inverted physical mass ordering scenarios and, in the case of a positive detection of neutrino masses, would recover the correct $\sum m_\nu$ value without bias [34, 164]. We also impose a minimal physical (uniform) prior $\sum m_\nu>0$ eV, while noting that scenarios with negative ‘effective’ neutrino mass have also recently been discussed (e.g., [161, 162, 163]). Supporting paper [50] explores the implications of different modelling assumptions in greater detail.

Primordial neutrinos were copiously produced before the nucleosynthesis era and decoupled at a temperature of about $1$ MeV; see, e.g., [34] for a review of neutrino cosmology. As the Universe expanded, neutrinos gradually lost kinetic energy, behaving as radiation in the early Universe and transitioning to non-relativistic matter around redshifts of $z\sim100$ for realistic neutrino masses, thereafter influencing the late-time expansion history by contributing to the matter component. The main effect of massive neutrinos on the CMB is to impact the angular diameter distance to last scattering, which is degenerate with the effects of other cosmological parameters such as $\Omega_m$ and $H_0$ (see, e.g., [165] for a recent discussion). Neutrinos also affect the lensing of CMB anisotropies by suppressing the growth of structure below the free-streaming scale. BAO are not sensitive to the latter effect at all, and only probe the background geometry by constraining the total matter density $\Omega_m$ and the parameter combination $H_0r_d$, so DESI BAO alone cannot constrain the neutrino masses. Nevertheless by breaking geometrical degeneracies, BAO significantly enhance the ability of the CMB to constrain this parameter.

The upper limits on $\sum m_\nu$ that we obtain from the combination of DESI and CMB depend on the particular choice of the CMB likelihood used, since the various likelihoods differ slightly in the amount of lensing power they infer from the lensed $TT$, $TE$ and $EE$ power spectra. This can be incorporated into a phenomenological parameter $A_L$ that scales the model lensing power used to compute the lensed power spectra (but not the power reconstructed from the 4-point function), such that values $A_L>1$ indicate an excess of lensing power, often referred to as the ‘$A_L$ anomaly’. Increasing $\sum m_\nu$ above $0.06$ eV decreases the expected lensing power, while smaller values would boost it.

The baseline CMB likelihood we have adopted through this paper combines the CamSpec likelihood [108, 109] together with the combined CMB likelihood from Planck [166] and ACT [104], both of which are based on the PR4 data release. Combining these with the DESI DR2

BAO, and assuming a $\Lambda$CDM background, we find the marginalized posterior limit

Although the CamSpec likelihood appears to prefer a lower lensing power and reduces the $A_L$ anomaly [109], the tighter constraints on other cosmological parameters degenerate with $\sum m_\nu$ mean that this upper bound is tighter than the one obtained with the original Plik likelihood:

The result of eq. (30) is directly comparable to the result obtained from BAO and CMB in the DESI DR1 analysis [38] [46], which used the Plik likelihood as the default and obtained the corresponding constraint $\sum m_\nu< 0.082$ eV (95%).²¹ This close to 20% reduction in the upper bound reflects the higher precision of the DESI DR2 BAO data compared to DR1.

The use of the alternative LoLLiPoP and HiLLiPoP (L-H) likelihoods of [110] also show a lower excess lensing power than Plik. In this case, this effect is also

reflected in significantly looser upper bounds on the neutrino mass,

(95%, DESI+CMB[L–H]).

Figure 15 shows the marginalized 1D posterior constraints on $\sum m_\nu$ for the combination of DESI with each of the CMB likelihoods.

Irrespective of which CMB likelihood is used, all three upper bounds are very close to the lower bound $\sum m_\nu> 0.059$ eV set by neutrino oscillation experiments [167]. The posteriors shown in Figure 15 all resemble the tails of distributions that would peak at negative neutrino masses if this were allowed by the prior. This behavior has been explored in a few recent studies [161, 163, 168, 169], as well as in our companion work [50], which also shows that it holds for profile likelihoods based on extrapolating the $\chi^2$ likelihood surface to $\sum m_\nu< 0$ eV. In the context of a model with effective neutrino masses that are allowed to assume negative values, the DESI+CMB likelihood peaks at $\sum m_{\nu,\mathrm{eff}} < 0$, and the tension with the minimal mass bound from terrestrial experiments increases to $3.0\sigma$ [50].

This discrepancy can be related to the tension between DESI and CMB results for the parameters $\Omega_m$ and $H_0r_d$ in a $\Lambda$CDM background discussed in Section VI above, and particularly in Figure 8. For fixed angular acoustic scale $\theta_*$, the effects on the CMB of changing $\sum m_\nu$ and changing $\Omega_m$ are positively correlated, as can be seen from the dashed black contour in Figure 16 (similar degeneracies can also be observed between $\sum m_\nu$ and $H_0r_d$, but with a negative correlation).

but with a negative correlation). In $\Lambda$CDM, DESI tightly constrains $\Omega_m$, as shown by the vertical shaded contour, to a value close to the low-end tail of the marginalized posterior distribution for CMB measurements. Consequently, the combination of DESI and CMB rules out all but the smallest neutrino masses. This result is thus related to the preference of DESI for lower values of $\Omega_m$ compared to the CMB discussed in the previous Sections.

As the upper limits on $\sum m_\nu$ from cosmology approach the lower limits from neutrino oscillations, the distinction between the different mass orderings becomes increasingly relevant. The baseline constraint eq. (29), obtained assuming a $\sum m_\nu> 0$ prior, already appears to rule out the inverted ordering. However, to determine the evidence in support of the normal ordering, one should also account for the fact that much of the posterior volume violates the constraints for both mass orderings. To do this we use a physical model for the total mass $\sum m_\nu$ in terms of lightest neutrino mass $m_\ell\ge0$, which includes information on the mass squared splittings $\Delta m^2_{21}$ and $\lvert\Delta m^2_{31}\rvert$ and allows for either mass ordering, as described in [170]. We then determine the upper bounds on this model from a combined analysis of DESI, CMB (with the baseline CamSpec likelihood) and 3-flavor oscillation constraints from NuFit-6.0 [167], assigning equal prior probabilities to either mass ordering. The result obtained is

(95%, DESI+CMB+NuFit-6.0).

The difference between the limits in eq. (29) and eq. (32) arises because the former only requires $\sum m_\nu> 0$ while the latter incorporates oscillation data that effectively require $\sum m_\nu> 0.06\,\mathrm{eV}$, therefore pushing the upper limit into territory disfavored by the DESI+CMB data alone. In the supporting paper [50], we also compare the two mass orderings and find a Bayes factor of $K = 10$ providing strong evidence in support of the normal mass ordering under the assumption of the $\Lambda$CDM+$\sum m_\nu$ model.

Conclusions

We have presented BAO measurements from over 14 million discrete galaxy and quasar tracers drawn from the first three years of operation of DESI and which will be included in the second data release (DR2). These results use samples of nearby bright galaxies, LRGs, ELGs and quasars over the redshift range $0.1 < z < 2.1$, and cover a cumulative effective volume of over 42 Gpc$^3$. Complementary BAO measurements from correlations of the Ly$\alpha$ forest and high-redshift quasars at effective redshift $z = 2.33$ are presented in the companion paper [63].

Our BAO analysis largely follows the methods used for the previous DESI DR1 analysis and presented in [21], but with improved statistical precision as the effective volume of the data has increased by more than a factor of two. Some particular differences include the use of a fainter limiting magnitude cut to define the bright galaxy sample, resulting in a higher number density, and the inclusion of quadrupole information in BAO fits to the quasar sample. As for the DR1 analysis, we applied a strict catalog-level blinding to our data while initial data checks were carried out and the analysis pipeline was being finalized. The validation tests and the criteria that were required to be met before the data were unblinded are described in detail in the supporting publication [48].

The final BAO results provide a precision on the isotropic distance scale measurement of $D_V/r_d$ that ranges from 1.54% (for the QSO) down to just 0.45% for our most constraining composite LRG3+ELG1 sample; other than for the QSO, the precision is sub-percent for every tracer and redshift range. We also obtained a precision of a few percent on the measurement of the Alcock-Paczyński distance ratio for every tracer except the lowest redshift BGS. Together these results are the most precise BAO measurements ever made at all redshifts covered, including at $z < 0.8$ where DESI DR2 now greatly exceeds the precision of the Sloan Digital Sky Survey (SDSS).

The DR2 results are very consistent with those previously reported from the DR1 data that form a subset of DR2. The $p$-value determined from the Kolmogorov-Smirnov test for the distribution of the statistics is 0.40. The DR2 results are also consistent with SDSS. [21] previously reported a $\sim3\sigma$ discrepancy in the measurement of the transverse BAO scale $D_M/r_d$ in the LRG2 sample at effective redshift $z_{\rm eff} = 0.706$ compared to the equivalent result from SDSS. For DR2, the difference has decreased in significance, lying within the range $1.5\sigma$ to $2.6\sigma$, depending on assumptions about the degree of correlation between the two samples.

The combination of all BAO measurements from DR2 is well fit by a flat $\Lambda$CDM cosmological model with matter density parameter $\Omega_m = 0.2975 \pm0.0086$ and product of the scaled Hubble constant and sound horizon at the drag epoch $h r_d = (101.54 \pm0.73)$ Mpc. This represents a 40% improvement in precision compared to the equivalent results from DR1, but with excellent consistency be-

Given the discrepancy between DESI and CMB under the $\Lambda$CDM model, and the strong preference for a $w_0w_a$CDM model discussed in Section VII, we consider the effect of allowing the equation-of-state parameters $w_0$ and $w_a$ to vary. In this scenario and assuming a $\sum m_\nu> 0$ prior, we find

from combining DESI with CMB, and

when also using DESY5 SNe to better constrain the equation of state of dark energy. More details of parameter posteriors in this scenario, including those obtained with other data combinations and the $w$CDM+$\sum m_\nu$ model, are presented in Table VII. Through comparison to Table V, we find that allowing $\sum m_\nu$ to vary with a conventional physical prior $\sum m_\nu> 0\,\mathrm{eV}$ has little effect on the $w(z)$ constraints, but that allowing $w_0$ and $w_a$ to vary significantly relaxes the neutrino mass bound in all cases. Moreover, in the case of DESI and CMB, the peak of the 1D marginalized posterior is recovered in the positive mass range, as shown in Figure 15.

Model/Dataset	$\Omega_m$	$H_0$ [km s$^{-1}$ Mpc$^{-1}$]	$H_0r_d$ [100 km s$^{-1}$]	$\sum m_\nu$ [eV]	$w$ or $w_0$	$w_a$
$\Lambda$CDM+$\sum m_\nu$
DESI BAO+CMB [Camspec]	0.3009 $\pm$ 0.0037	68.36 $\pm$ 0.29	100.96 $\pm$ 0.48	$< 0.0642$	—	—
DESI BAO+CMB [L-H]	0.2995 $\pm$ 0.0037	68.48 $\pm$ 0.30	101.16 $\pm$ 0.49	$< 0.0774$	—	—
DESI BAO+CMB [Plik]	0.2998 $\pm$ 0.0038	68.56 $\pm$ 0.31	101.09 $\pm$ 0.50	$< 0.0691$	—	—
$w$CDM+$\sum m_\nu$
DESI BAO+CMB	0.2943 $\pm$ 0.0073	69.28 $\pm$ 0.92	102.3 $\pm$ 1.3	$< 0.0851$	$-1.039 \pm0.037$	—
DESI BAO+CMB+Pantheon+	0.3045 $\pm$ 0.0051	67.94 $\pm$ 0.58	100.35 $\pm$ 0.84	$< 0.0653$	$-0.985 \pm0.023$	—
DESI BAO+CMB+Union3	0.3047 $\pm$ 0.0059	67.93 $\pm$ 0.69	100.33 $\pm$ 0.99	$< 0.0649$	$-0.985 \pm0.028$	—
DESI BAO+CMB+DESY5	0.3094 $\pm$ 0.0049	67.34 $\pm$ 0.53	99.49 $\pm$ 0.78	$< 0.0586$	$-0.961 \pm0.021$	—
$w_0w_a$CDM+$\sum m_\nu$
DESI BAO+CMB	0.353 $\pm$ 0.022	$63.7^{+1.7}_{-2.2}$	$93.8^{+2.5}_{-3.2}$	$< 0.163$	$-0.42^{+0.24}_{-0.21}$	$-1.75 \pm0.63$
DESI BAO+CMB+Pantheon+	0.3109 $\pm$ 0.0057	67.54 $\pm$ 0.59	99.62 $\pm$ 0.86	$< 0.117$	$-0.845 \pm0.055$	$-0.57^{+0.23}_{-0.19}$
DESI BAO+CMB+Union3	0.3269 $\pm$ 0.0088	65.96 $\pm$ 0.84	97.3 $\pm$ 1.2	$< 0.139$	$-0.674 \pm$ 0.090	$-1.06^{+0.34}_{-0.28}$
DESI BAO+CMB+DESY5	0.3188 $\pm$ 0.0058	66.75 $\pm$ 0.56	98.43 $\pm$ 0.83	$< 0.129$	$-0.758 \pm$ 0.058	$-0.82^{+0.26}_{-0.21}$

In $\Lambda$CDM, BAO measurements from DESI set a geometric constraint on $\Omega_m$ (or, equivalently, on $H_0r_d$) that helps break the degeneracy in the CMB results for $\sum m_\nu$. In the $w_0w_a$CDM model, these geometric measurements primarily help by constraining the background $w_0$ and $w_a$ values, thus limiting the available $\sum m_\nu$ range from the CMB alone, which is itself weaker than in $\Lambda$CDM due to the additional freedom in the expansion history. This effect is even clearer when also adding SNe to DESI and CMB, as the constraints on the background become significantly tighter. This is illustrated in Figure 17, where the vertical blue bands show the purely geometrical constraints from DESI+DESY5. An extension of this plot where these degeneracies extend further in parameter space by allowing for negative neutrino masses through and effective parameter $\sum m_{\nu,\mathrm{eff}} > 0$ is shown in [50].

Finally, constraints on the number of effective relativistic degrees of freedom, $N_{\rm eff}$, are not discussed in detail here, but are consistent with the particle physics value $N_{\rm eff} = 3.044$ as described in our supporting paper [50]. between the two. However, while a $\Lambda$CDM model provides a good fit to DESI data, the model parameters obtained from this fit are now in $2.3\sigma$ tension with those derived from the CMB, increased from $1.9\sigma$ in DR1. This tension is present despite DESI being consistent with the acoustic angular scale $\theta_\ast$ measured by CMB.

When calibrated with an external prior on $\Omega_{\mathrm{b}}h^2$ from BBN, our BAO measurements correspond to a value $H_0 = (68.50 \pm0.58)\ \mathrm{km\,s^{-1}\,Mpc^{-1}}$ in $\Lambda$CDM, a value that is independent of any information on CMB anisotropies. In the $\Omega_{\mathrm{m}}$-$H_0$ plane, the results from BBN-calibrated BAO are now more discrepant with those from CMB; the offset of the results is again along the $\Omega_{\mathrm{m}}h^3 = \mathrm{constant}$ degeneracy direction of the CMB that is determined by the very precisely measured acoustic angular scale $\theta_\ast$. Combining DESI BAO with BBN and the Planck $\theta_\ast$ result gives $H_0 = (68.45 \pm0.47)\ \mathrm{km\,s^{-1}\,Mpc^{-1}}$, a 0.7% precision measurement competitive with that from the CMB itself, and in strong disagreement with SH0ES [129].

Within the $\Lambda$CDM framework, DESI results are also somewhat in tension with the high $\Omega_{\mathrm{m}}$ values preferred by SNe datasets, which—contrary to DESI—prefer larger $\Omega_{\mathrm{m}}$ than Planck. While not individually rising to the $3\sigma$ significance threshold, these results point to an incompatibility between different cosmological datasets when interpreted in the $\Lambda$CDM model. Interestingly, the relative levels of $\Omega_{\mathrm{m}}$ and $H_0$ values currently measured by SN, DESI and CMB datasets match what would be expected if data from a true evolving dark energy model were analyzed in the restrictive $\Lambda$CDM model, as pointed out recently by [131].

Assuming a $\Lambda$CDM background, the combination of DESI and CMB data give the tightest upper bound on the neutrino mass sum to date, $\sum m_\nu< 0.064\ \mathrm{eV}$ (95% limit) in our baseline analysis. Although this relaxes to $\sum m_\nu< 0.078\ \mathrm{eV}$ when using an alternative Planck likelihood [110], both results are approaching the lower bound set by terrestrial neutrino oscillation experiments, $\sum m_\nu\ge0.059\ \mathrm{eV}$. Indeed, if the model is extended to allow for negative ‘effective’ neutrino masses, most of the posterior mass and the peak of the likelihood lie at $\sum m_{\nu,\mathrm{eff}} < 0\ \mathrm{eV}$, another possible sign of growing tensions within the $\Lambda$CDM model with DESI DR2 [50].

The evidence for a departure from $\Lambda$CDM in the form of evolving dark energy has increased with the DR2 BAO data. Comparing the evolving dark energy model parametrized by $w(a) = w_0 + w_a(1-a)$ to $\Lambda$CDM, we find a $3.1\sigma$ evidence in favor of dynamical dark energy from DESI+CMB alone. When we add the recent Pantheon+, Union3 or DESY5 SNe datasets to this combination, the preference for $w_0w_a$CDM over $\Lambda$CDM is $2.8\sigma$, $3.8\sigma$ or $4.2\sigma$ respectively, with all three giving results for $w_0$ and $w_a$ consistent with each other within their 68% credible regions. The preference for $w_0w_a$CDM under these dataset combinations have increased compared to our DR1 results in [38]. All combinations of available datasets we study favor mutually consistent results with $w_0 > -1$ and $w_a < 0$, apparently indicating a weakening dark energy today and a phantom crossing at some point in the past (in this parametrization). The supporting paper [49] examines this behavior using a wider range of models and non-parametric methods. We note that the degeneracy direction for the constraints in the $w_0$-$w_a$ plane approximately points towards the $\Lambda$CDM solution, although this is neither exact nor consistent between fits to different combinations of datasets. This constitutes a weak coincidence as there is no a priori reason for the best-fit values of $w_0$ and $w_a$ to have the values that lead to this observation.

In Sections VIII B and VIII C, we examined the contributions to the preference for evolving dark energy in detail. For DESI+CMB, the discrepancy in preferred $\Omega_{\mathrm{m}}$ values in $\Lambda$CDM is resolved in $w_0w_a$CDM, leading to the $3.1\sigma$ ($\Delta\chi^2_{\mathrm{MAP}} = -12.5$) preference for the latter. If we omit CMB lensing, the best-fit $w_0w_a$CDM parameters are unchanged, though the significance drops to $2.7\sigma$ ($\Delta\chi^2_{\mathrm{MAP}} = -9.7$). If we also marginalize over any late-time dependence in the CMB and only impose a prior on $(\theta_\ast,\omega_{\mathrm{b}},\omega_{\mathrm{c}})$, there is a further slight reduction to $2.4\sigma$ ($\Delta\chi^2_{\mathrm{MAP}} = -8.0$). The $w_0w_a$CDM and ‘compromise’ $\Lambda$CDM models predict different evolution of the distance-redshift relation at $z < 0.4$, with a $\sim2\%$ difference in distances at $z \sim0.1$. BAO do not provide high precision in this regime, but supernovae do, and the DESY5 and Union3 samples clearly favor the $w_0w_a$CDM prediction. However, the distance-redshift evolution in the Pantheon+ analysis is closer to the $\Lambda$CDM prediction, so adding Pantheon+ to DESI+CMB does not strengthen the evidence for $w_0w_a$CDM. Appendix A shows that our results are robust to changes in the choice of CMB likelihood. We also show that replacing the CMB likelihood with the weak gravitational lensing constraints from DESY3 (3$\times$2pt) analysis retains a clear preference for evolving dark energy, with, e.g., a $\gtrsim3\sigma$ preference for $w_0w_a$CDM from DESI+DESY3 (3$\times$2pt)+DESY5 SNe. Removing the preference for $w_0w_a$CDM through changes to the DESI measurements themselves would require systematic errors to be far larger than any found in our tests.

In the $w_0w_a$CDM background model, when the neutrino mass scale is allowed to vary with physical non-negative prior bounds, the upper limit on $\sum m_\nu$ is significantly relaxed to $\sum m_\nu< 0.16\ \mathrm{eV}$ (95%, from DESI+CMB) or $\sum m_\nu< 0.13\ \mathrm{eV}$ (95%, from DESI+CMB+DESY5), the middle result from the three possible SNe samples), while the constraints on $w_0$ and $w_a$ do not materially change from the case where $\sum m_\nu$ is fixed. These limits are entirely consistent with neutrino oscillation experiment results. The supporting paper [50] shows that even when allowing for effective masses $\sum m_{\nu,\mathrm{eff}}$ with a wide prior that allows negative values, the marginalized posterior in $w_0w_a$CDM from DESI+CMB is consistent with the oscillation lower limit and its peak lies in the positive region.

In [38] we characterized the results from DESI DR1 and external datasets as providing “a tantalizing suggestion of deviations from the standard cosmological model”. The DR2 data presented here have sharpened this evidence, although significance levels still vary depending on the external data used—particularly the choice of SNe sample—and no combination exceeds $4.2\sigma$. Nevertheless, it is becoming clear that unless some unidentified systematic error affects one or several of the different cosmological datasets used, the challenge to the $\Lambda$CDM model has increased. Sharper measurements from future DESI analyses and from other experiments will show whether these challenges to the standard cosmological model herald yet another radical transformation in our understanding of the evolution and energy content of the Universe.

Data Availability

The data used in this analysis will be made public along with Data Release 2 (details in https://data.desi.lbl.gov/doc/releases/). BAO likelihoods for DESI DR2, integrated in the cobaya code, are available at https://github.com/CobayaSampler/bao_data. The data points corresponding to the figures from this paper are available in Ref. [171].

Acknowledgments

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 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.

The DESI Legacy Imaging Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS), and the Mayall z-band Legacy Survey (MzLS). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory. Legacy Surveys also uses data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Legacy Surveys was supported by: the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility; the U.S. National Science Foundation, Division of Astronomical Sciences; the National Astronomical Observatories of China, the Chinese Academy of Sciences and the Chinese National Natural Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. The complete acknowledgments can be found at https://www.legacysurvey.org/.

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 I’oligam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

References

Early-Universe priors and robustness to choice of CMB likelihoods

As discussed in Section IV, for some purposes it is convenient to make use of a compression of the full CMB information into a multivariate correlated Gaussian prior on the quantities $\theta_*$, $\omega_{\rm bc}$ and $\omega_b$, which are early-Universe quantities that can be determined independently of assumptions about the late-time evolution by marginalizing over late-time effects such as the ISW effect and CMB lensing. This process was performed by [111] based on the CamSpec CMB likelihood. We use their results to define a Gaussian prior with mean

and covariance

The use of this prior has been indicated throughout the paper by the use of the $(\theta_*,\omega_b,\omega_{\rm bc})_{\rm CMB}$ notation. As shown in the panel of Figure 14 the posterior results for the $w_0w_a$CDM model are remarkably similar in all cases when comparing fits to DESI+$(\theta_*,\omega_b,\omega_{\rm bc})_{\rm CMB}$, DESI+CMB (with full CMB information) and DESI+CMB (no lensing), indicating that most of the information added by the CMB in this scenario comes from the anchoring constraint it provides on $\Omega_m$, limiting the freedom of the model to absorb background expansion effects of $w_0$ and $w_a$ into other parameters. We also tested a separate compression of the CMB information based on the shift parameters $R$ and $l_A$ together with $\omega_b$ [112] and found that it gave very consistent results as well (see also [172]). Nevertheless, while the addition of CMB lensing information does not noticeably shift the $w_0$-$w_a$ posteriors, it does appreciably change the $\Delta\chi^2_{\rm MAP}$ values and thus the significance of rejection of $\Lambda$CDM.

Figure 18 shows the differences in the constraints obtained from the DESI+CMB+DESY5 data combination in the $w_0w_a$CDM case, arising from the use of the three different CMB likelihoods, namely CamSpec (our baseline in this paper), Plik and L–H (cf. Section IV for descriptions). The results using CamSpec and Plik are remarkably consistent, though we see a small shift at the $0.1\sigma$ level when using L–H PR4.

Robustness to assumption about correlation in BAO systematic errors

Our baseline approach, described briefly in § 3 and in more detail in [48], assumes that the systematic error contributions to BAO measurements are not correlated with each other across redshift bins. As a robustness test, we evaluated the impact of introducing off-diagonal terms in the systematics covariance matrix of the BAO measurements. To do this, we assumed that the theoretical modelling systematic error contributions, which are estimated to be 0.1% in $\alpha_{\rm iso}$ and 0.2% in $\alpha_{\rm AP}$, are correlated across all redshift bins, with a correlation coefficient of 0.5. This leads to redshift-dependent correlations in the final covariance matrix. We found that including these correlations had no detectable effect on the cosmological constraints obtained in any model or combination of DESI with any external dataset. This result is unsurprising, since the systematic errors are subdominant to statistical errors in all redshift bins, so correlations in the systematic errors cannot play a large role.