This paper provides a comprehensive overview of how fitting of baryon acoustic oscillations BAO is carried out within the upcoming Dark Energy Spectroscopic Instruments DESI 2024 results using its DR1 data set and the associated systematic error budget from theory and modelling of the BAO. We derive new results showing how nonlinearities in the clustering of galaxies can cause potential biases in measurements of the isotropic inlineformulatexmath idTM0001 notationLaTeXalpha mathrmisotexmathinlineformula and anisotropic inlineformulatexmath idTM0002 notationLaTeXalpha mathrmaptexmathinlineformula BAO distance scales and how these can be effectively removed with an appropriate choice of reconstruction algorithm. We then demonstrate how theory leads to a clear choice for how to model the BAO and develop implement and validate a new model for the remaining smoothbroadband i.e. without BAO component of the galaxy clustering. Finally we explore the impact of all remaining modelling choices on the BAO constraints from DESI using a suite of highprecision simulations arriving at a set of best practices for DESI BAO fits and an associated theory and modelling systematic error. Overall our results demonstrate the remarkable robustness of the BAO to all our modelling choices and motivate a combined theory and modelling systematic error contribution to the postreconstruction DESI BAO measurements of no more than 0.1 per cent 0.2 per cent for its isotropic anisotropic distance measurements. We expect the theory and best practices laid out to here to be applicable to other BAO experiments in the era of DESI and beyond.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Baryon acoustic oscillation theory and modelling systematics for the DESI 2024 results
{'Baryon Acoustic Oscillation Theory and Modelling Systematics for the DESI 2024 results': 'Shi-Fan Chen 1 ★ ‡ , Cullan Howlett 2 †‡ , Martin White, 3 , 4 Patrick McDonald, 4 Ashley J. Ross, 5 Hee-Jong Seo, 6 Nikhil Padmanabhan, 7 J. Aguilar, 4 S. Ahlen, 8 S. Alam, 9 O. Alves, 10 U. Andrade, 11 , 10 R. Blum, 12 D. Brooks, 13 X. Chen, 7 S. Cole, 14 K. Dawson, 15 A. de la Macorra, 16 Arjun Dey, 12 Z. Ding, 17 P. Doel, 13 S. Ferraro, 4 , 18 A. Font-Ribera, 13 , 19 D. Forero-Sánchez, 20 J. E. Forero-Romero, 21 , 22 C. Garcia-Quintero, 23 , 24 E. Gaztañaga, 25 , 26 , 27 S. Gontcho A Gontcho, 4 M. M. S Hanif, 10 K. Honscheid, 5 , 28 , 29 T. Kisner, 4 A. Kremin, 4 A. Lambert, 4 M. Landriau, 4 M. E. Levi, 4 M. Manera, 30 , 19 A. Meisner, 12 J. Mena-Fernández, 31 R. Miquel, 32 , 19 A. Munoz-Gutierrez, 16 E. Paillas, 33 , 34 N. Palanque-Delabrouille, 35 , 4 W. J. Percival, 33 , 36 , 34 A. Pérez-Fernández, 16 , 37 F. Prada, 38 M. Rashkovetskyi, 23 M. Rezaie, 39 A. Rosado-Marin, 6 G. Rossi, 40 R. Ruggeri, 41 , 2 E. Sanchez, 42 D. Schlegel, 4 J. Silber, 4 G. Tarlé, 10 M. Vargas-Magaña, 16 B. A. Weaver, 12 J. Yu, 20 S. Yuan, 43 R. Zhou, 4 and Z. Zhou 44  \nAffiliations listed in Appendix G.  \nAccepted 2024 September 3. Received 2024 September 3; in original form 2024 February 22', 'ABSTRACT': "This paper provides a comprehensive overview of how fitting of Baryon Acoustic Oscillations (BAO) is carried out within the upcoming Dark Energy Spectroscopic Instrument's (DESI) 2024 results using its DR1 dataset, and the associated systematic error budget from theory and modelling of the BAO. We derive new results showing how non-linearities in the clustering of galaxies can cause potential biases in measurements of the isotropic ( 𝛼 iso ) and anisotropic ( 𝛼 ap ) BAO distance scales, and how these can be effectively removed with an appropriate choice of reconstruction algorithm. We then demonstrate how theory leads to a clear choice for how to model the BAO and develop, implement and validate a new model for the remaining smooth-broadband (i.e., without BAO) component of the galaxy clustering. Finally, we explore the impact of all remaining modelling choices on the BAO constraints from DESI using a suite of high-precision simulations, arriving at a set of best-practices for DESI BAO fits, and an associated theory and modelling systematic error. Overall, our results demonstrate the remarkable robustness of the BAO to all our modelling choices and motivate a combined theory and modelling systematic error contribution to the post-reconstruction DESI BAO measurements of no more than 0 . 1% (0 . 2%) for its isotropic (anisotropic) distance measurements. We expect the theory and best-practices laid out to here to be applicable to other BAO experiments in the era of DESI and beyond.  \nKey words: cosmology: theory, large-scale structure of Universe, distance scale, cosmological parameters", '1 INTRODUCTION': 'Constraining the expansion history of the Universe is important both for improving our empirical cosmographical knowledge and for the constraints it places on the constituents of the Universe and our theories of gravity. Within the standard Friedmann-Lemaître-Robertson-Walker (FLRW; Friedmann 1922; Lemaître 1931; Robertson 1935; Walker 1937) model, the scale factor 𝑎 ( 𝑡 ) - is the only degree of freedom at the background level, determining distance relations over cosmic history. Constraining its behaviour has been a major goal of cosmology for close to a century.  \n- ★ E-mail: sfschen@ias.edu\n- † E-mail: c.howlett@uq.edu.au\n- ‡ Both authors contributed equally to this work.  \nOne of the most robust and accurate ways of determining 𝑎 ( 𝑡 ) is through measurements of baryon acoustic oscillations (BAO) in the distribution of galaxies and gas. Oscillations in the baryon-photon fluid prior to decoupling leave an imprint in both the cosmic microwave background and the matter density with a characteristic length scale, related to the distance 𝑟 𝑑 a sound wave can travel prior to the baryon-drag epoch, just after recombination. Measurements of this scale at low redshift in large-scale structure provide a \'standard ruler\' for constraining the expansion history of the Universe, both by measuring the angular diameter distance ( 𝐷 𝐴 ( 𝑧 ) ) and the Hubble parameter ( 𝐻 ( 𝑧 ) ) to a specific redshift. This ruler can be made \'absolute\' by connecting the BAO scale back to early universe physics, e.g. through constraints from the cosmic microwave background (CMB; Peebles & Yu 1970; Sunyaev & Zeldovich 1970; Dodelson & Schmidt  \n2020; Huterer 2023) or big-bang nucleosynthesis (BBN; Cooke et al. 2018).  \nThe BAO method has become a workhorse of modern cosmology and provides some of our tightest current constraints on the cosmological expansion history (Cole et al. 2005; Eisenstein et al. 2005; Beutler et al. 2011; Blake et al. 2011; Alam et al. 2021). Relying as it does on a feature at large spatial separations that is little affected by non-linear evolution or the complex astrophysics of galaxy formation, it is robust and theoretically well understood. For this reason, it is generally regarded as a low-systematics method for constraining cosmology. Indeed, significant theoretical effort over the past two decades has led to a fairly mature understanding of the effect of the small shifts and modulations to the BAO signal incurred by largescale structure physics which we recap and build on in this work (§2).  \nThe purpose of this paper is to use our theoretical understanding to generate a theory and modelling systematics budget, and suggest fitting best practices, for BAO measurements in the Dark Energy Spectroscopic Instrument (DESI). DESI is a Stage-IV spectroscopic instrument (DESI Collaboration et al. 2016, 2022) that will deliver galaxy BAO measurements with unprecedented precision up to 𝑧 ≈ 2 from galaxy clustering alone, with even greater reach when Lyman𝛼 data are included at higher redshift. At the end of its five-year run DESI is expected to yield measurements on the (isotropic) BAO scale at a cumulative 0 . 2% precision (DESI Collaboration et al. 2016, 2024d). 1 Already, galaxies from the early data release (EDR) in the luminous red galaxy (LRG) sample and bright galaxy sample (BGS), representing a small fraction of the total DESI data, yield BAO detections at the few percent level (DESI Collaboration et al. 2024e; Moon et al. 2023).  \nThis paper is part of a set focused on the analysis of clustering in DESI Data Release 1 (DR1; DESI Collaboration 2024a,b), and in support of the key paper presenting the main galaxy BAO measurements for that data set (DESI Collaboration et al. 2024a). The DR1 data further contains BAO information in the Lyman𝛼 forest (DESI Collaboration et al. 2024b), and the cosmological implications of the joint galaxy and Lyman𝛼 BAO measurements are discussed in DESI Collaboration et al. (2024c). Information beyond the BAO - particularly the fullshape information in the galaxy 2-point function-are studied in DESI Collaboration (2024c,d,e). While the aggregate isotropic BAO scale precision from analysis of DR1 is closer to 0 . 5%, our conclusions and recommendations are somewhat general and aimed towards the full DESI sample. Our understanding of galaxy clustering in the standard model of cosmology ( Λ CDM) gives us a baseline of known theoretical systematics that need to be considered in order to make robust measurements of the BAO scale, even before unknown modifications due to nonstandard physics are considered.  \nThe structure of the paper is as follows: In Section 2 we begin with a broad overview of the BAO signal as observed in large-scale structure simulations. In Sections 3 and 4 we review the theory of the BAO and BAO reconstruction, computing error estimates due to theory systematics, filling in extant gaps in the literature where necessary, and motivating our baseline form for fitting the BAO signal. In Section 5 we introduce a cubic-spline based method for  \n1 Where the cumulative or "aggregate" BAO precision is the statistical uncertainty we would expect in the BAO scale if all DESI galaxies were put into the same redshift bin and used to produce a single BAO measurement. It is equivalent to an inverse sum of the BAO precision across all our independent measurement bins.  \nmarginalizing over the non-BAO broadband of galaxy clustering and harmonize broadband models in Fourier and configuration space. Section 6 introduces the N-body simulations and inference pipeline B/a.pc/r.pc/r.pc/y.pc implementing the suggested analysis choices in this paper, and the systematic effects of various numerical choices are examined in Section 7. We conclude in Section 8.', '2 THE OBSERVED BARYON ACOUSTIC OSCILLATION SIGNAL: AN OVERVIEW': 'Baryon acoustic oscillations (BAO) manifest as an oscillatory, or wiggly, imprint in the linear matter power spectrum 𝑃 lin ( 𝑘 ) . The physics is very well understood and can be computed to high accuracy by numerical evolution of the coupled Einstein, fluid and Boltzmann equations as is done e.g. by CAMB (Lewis et al. 2000; Howlett et al. 2012) and CLASS (Blas et al. 2011). To recap the basic physics, prior to recombination, when the densities are high, baryon-photon scattering is rapid compared with the travel time across a wavelength. This allows the behaviour of the baryon and photons to be described as a single fluid, whose continuity and Euler equations reduce to a driven harmonic oscillator with slowly varying natural frequency. During this \'tight-coupling\' phase the baryon perturbations hence also undergo harmonic motion with very slowly decaying amplitude (Peebles & Yu 1970; Sunyaev & Zeldovich 1970). These \'baryon acoustic oscillations\' become \'frozen-in\' at decoupling when the photons release their hold on the baryons, setting the oscillation frequency 𝑟 𝑑 , and the matter perturbations begin to grow as ∝ 𝑎 . The combination of oscillatory densities and velocities which project onto the growing mode gives the final post-recombination spectrum. 2 This then grows in a scale-independent fashion until the present day, up to small corrections due to massive neutrinos which start off relativistic but transition to non-relativistic at some point post-recombination. The baryon acoustic oscillations themselves are superposed upon the smooth power spectrum, whose shape is primarily determined by the epoch of matter-radiation equality. The amplitude of the oscillations depends upon the driving force (gravitational potentials) and the baryon-to-photon ratio. Larger oscillations come from higher baryon density 𝜔 𝑏 and/or smaller total-matter density 𝜔 𝑚 . Additionally, light species (e.g. neutrinos) which propagate ahead of the sound waves also exert a gravitational influence on the BAO which can shift their phase (Bashinsky & Seljak 2004), and any isocurvature modes in the initial conditions can also alter the shape of the observed signal (Zunckel et al. 2011); the impact of some these effects on DESI BAO measurements is discussed in Pérez-Fernández et al. (2024), but we note that they are also more directly constrained by the CMB itself.  \nGiven a linear matter power spectrum 𝑃 lin , the redshift-space galaxy power spectrum in linear theory and in the plane-parallel or distant-observer approximation, is given by (Kaiser 1987)  \n𝑃 𝑔𝑔 ( 𝑘, 𝜇 ) = ( 𝑏 1 + 𝑓 𝜇 2 ) 2 𝑃 lin ( 𝑘 ) . (1)  \nwhere 𝑓 = 𝑑 ln 𝐷 / 𝑑 ln 𝑎 ≈ Ω 0 . 55 𝑚 is the growth rate and 𝑏 1 is the linear (scale-independent, deterministic) bias. Since the linear spectrum is otherwise smooth, and neither the nonlinearities of structure/galaxy formation nor survey systematics are expected to produce features with the same frequency, the aim of BAO measurements in  \n2 At low 𝑘 this is a mixture of density and velocity modes but at higher 𝑘 it is primarily the velocity mode which is why the higher harmonics of the BAO are out of phase with the peaks in the CMB, the latter of which arise primarily from photon density perturbations.  \ngalaxy surveys is to independently extract the oscillatory part of this linear signal from the overall clustering of galaxies and use it as a standard ruler. The BAO feature is a localized peak in the correlation function at a scale determined by the sound horizon at the baryondrag epoch ( 𝑟 𝑑 ), or equivalently as a series of oscillations in the power spectrum with Fourier-frequency 2 𝜋 / 𝑟 𝑑 . At a given redshift, 𝑧 , the line-of-sight and transverse determinations of this ruler fix 𝐻 ( 𝑧 ) 𝑟 𝑑 and 𝑟 𝑑 / 𝐷 𝐴 ( 𝑧 ) , respectively. Modern spectroscopic surveys typically convert the measured redshifts and sky positions of galaxies into comoving coordinates using a fiducial cosmology (usually flat Λ CDM with a fiducial Ω 𝑚, fid ≈ 0 . 3), in which the multipoles of the power spectrum 𝑃 ( 𝑘, 𝜇 ) = " ℓ 𝑃 ℓ ( 𝑘 )L ℓ ( 𝜇 ) and correlation functions 𝜉 ( 𝑟, 𝜇 ) can be computed in physical units. These 2-point functions are then fitted against a BAO \'template\' that is scaled in the line-of-sight and transverse directions (Padmanabhan & White 2008) to extract the so-called dilation parameters  \n𝛼 | | = 𝐻 fid ( 𝑧 ) 𝑟 tem 𝑠 𝐻 ( 𝑧 ) 𝑟 𝑠 , 𝛼 ⊥ = 𝐷 𝐴 ( 𝑧 ) 𝑟 tem 𝑠 𝐷 fid 𝐴 ( 𝑧 ) 𝑟 𝑠 , (2)  \nwhere \'fid\' and \'tem\' refer to quantities in the fiducial and template cosmologies respectively. It is also common to see the BAO dilation parameters recast into a purely isotropic term 𝛼 iso ≡ 𝛼 1 / 3 | | 𝛼 2 / 3 ⊥ and an anisotropic term 𝛼 ap ≡ ( 1 + 𝜖 ) 3 ≡ 𝛼 | | / 𝛼 ⊥ . We will use these all interchangeably in this work.  \nLet us describe the observation and extraction of the BAO in galaxy surveys outlined above in more detail, beginning with the dynamical model. In order to do so it is helpful to split the power spectrum into a BAOtemplate consisting of the wiggles alone 𝑃 𝑤 and a \'no-wiggle\' portion 𝑃 𝑛𝑤 containing the broadband power which we shall discuss further below (§7.4)  \n𝑃 lin ( 𝑘 ) = 𝑃 𝑛𝑤 ( 𝑘 ) + 𝑓 𝑏 𝑃 𝑤 ( 𝑘 ) . (3)  \nIn the definition of the wiggle/no-wiggle split above we have included in the wiggly component a pre-factor of the baryon fraction 𝑓 𝑏 ≈ 15%.Thissmallparametercontrolsthe size of the BAO signal relative to the overall clustering signal, and we will use it as a convenient power-counting parameter to denote the suppression of contributions involving higher powers of the wiggle component. However, since 𝑓 𝑏 is degenerate with the amplitude of 𝑃 𝑤 we will generally include it in 𝑃 𝑤 for convenience in the rest of the paper, except when it is useful for book-keeping like here. The measurement of the linear BAO signal can then be thought of as the measurement of the linear dependence of the observed galaxy power spectrum, 𝑃 𝑔𝑔 , on 𝑓 𝑏 . 3 Schematically, for the observed 2-point function we can write  \n𝑃 𝑔𝑔 ( 𝑘, 𝜇 ) = ( 𝑃 𝑔𝑔 ) 𝑓 𝑏 = 0 + 𝑓 𝑏 GLYPH<18> 𝑑𝑃 𝑔𝑔 𝑑𝑓 𝑏 GLYPH<19> 𝑓 𝑏 = 0 + O( 𝑓 2 𝑏 ) . (4)  \nBased on this structure models of the BAO typically take the form  \n𝑃 𝑔𝑔 ( 𝑘, 𝜇 ) = B( 𝑘, 𝜇 ) 𝑃 𝑛𝑤 ( 𝑘 ) + C( 𝑘, 𝜇 ) 𝑃 𝑤 ( 𝑘 ) + D( 𝑘, 𝜇 )  \nwhere the term proportional to C( 𝑘, 𝜇 ) nal 𝑑𝑃 𝑔𝑔 / 𝑑𝑓 𝑏 = 𝐶 ( 𝑘, 𝜇 ) 𝑃 𝑤 we want to measure, with 𝐶 ( 𝑘, 𝜇 )  \n(5) is the leading BAO sig-=  \n( 𝑏 + 𝑓 𝜇 2 ) 2 , in linear theory. The rough non-oscillatory broadband shape is contained in the B( 𝑘, 𝜇 ) term, given by the Kaiser formula in linear theory. The final term D( 𝑘, 𝜇 ) represents any remaining contributions to the observed power spectrum, which are largely due to nonlinear coupling of the smooth linear spectrum and observational systematics. It may also absorb higher order (e.g. ( 𝑓 𝑏 ) 2 ) terms in the BAO, though these will tend to be both rather suppressed and oscillatory. The systematics we will discuss in this paper largely have to do with the goodness of this approximation; specifically we will consider the extent to which nonlinear oscillatory terms modify the C( 𝑘, 𝜇 ) signal we wish to measure, and how well the remaining terms can be absorbed by D( 𝑘, 𝜇 ) under the assumption of smoothness.  \nGiven one adopts a fiducial cosmology for converting galaxy redshifts to distances, and a template linear spectrum with corresponding 𝑟 tem 𝑑 and 𝑃 tem 𝑤 , the assumption is then that any of the scale information carried in the BAO wiggles can be accounted for with a simple scaling 𝑃 𝑤 ( 𝑘 ) ∼ 𝑃 tem 𝑤 ( 𝑟 𝑑 𝑘 / 𝑟 tem 𝑑 ) in Equation 5. Allowing for 𝐻 ( 𝑧 ) and 𝐷 𝐴 ( 𝑧 ) to vary independently away from the fiducial cosmology, such that  \nk obs = 𝐷 𝐴 ( 𝑧 ) 𝐷 fid 𝐴 ( 𝑧 ) k true ⊥ , 𝐻 fid ( 𝑧 ) 𝐻 ( 𝑧 ) 𝑘 true ∥ ! (6)  \nleads to the conclusion that the free parameters required to \'match\' the template to the observed BAO are precisely the dilation parameters defined in Equation 2, with smooth changes to the broadband and the shape and normalization of the wiggles captured by marginalizing over B , C , D . This borrows philosophically from the AlcockPaczynski effect (Alcock & Paczynski 1979), which describes the dilation of scales and angles when analysing galaxy positions assuming a fiducial cosmology that differs from the true underlying cosmology of the Universe.  \nIt is worth noting that not all of the cosmological dependence of the BAO signal (e.g. its amplitude, damping or phase shifts) is captured by this rescaling - indeed, some of these, such as the phase shift,cannot even be captured by other parameters like galaxy bias or broadband modelling in standard BAO fits. In general, the physical effects that change 𝑃 lin also change the CMB anisotropies, often more significantly (Eisenstein & White 2004), and tend to be disfavored by existing measurements (e.g. Planck Collaboration et al. 2020). The systematic errors incurred by these additional cosmologydependent effects are characterized in an accompanying DESI paper (Pérez-Fernández et al. 2024), so we will only note here that the extraction of 𝛼 ∥ , ⊥ in BAO fits has been shown to be rather robust to these (Vargas-Magaña et al. 2014; Bernal et al. 2020; Carter et al. 2020), while referring interested readers to that upcoming companion paper for further details in the context of DESI Y1. In addition, the definition of 𝑃 𝑤 is not unique and several different template extraction methods exist in the literature. Properly constructed, these 𝑃 𝑤 should only differ by smooth functions of 𝑘 that are degenerate with the marginalized broadband terms, and we will explore the extent to which this is true in §7.4.  \nOur interest in the rest of this paper will be in quantifying the effects of dynamical nonlinearities and modelling choices on measuring 𝛼 ∥ , ⊥ so as to form a theoretical error budget for DESI BAO analyses. In doing so we will work primarily in the plane-parallel approximation and validate our models against the DESI Y1 \'cubic mocks\' based on single-redshift snapshots of N-body simulations in periodic, cubic boxes. This isolates the \'theoretical systematics\' from those depending upon observing geometry or observational non-idealities. However, before moving into these investigations it is important to describe a few additional details necessary to make contact with observations.  \nFirstly, galaxy surveys measure clustering not in translationally invariant cubic boxes but within defined angular and redshift windows. In the power spectrum this effect, along with wide-angle contributions to common power spectrum estimators, can be accounted for by convolving theory predictions with suitably computed window matrices (Feldman et al. 1994; Blake et al. 2018; Beutler & McDonald 2021). Since these effects should ideally be exactly computable given the galaxy window function we will not consider their systematics here but relegate them to the DESI Y1 catalog papers. Secondly, the measured 2-point function is not sampled from a single point but rather a weighted sum across the redshift range 𝑃 = " 𝑖 𝑤 𝑖 𝑃 ( 𝑧 𝑖 ) ; to leading order this is taken into account by taking the BAO measurement to be at the effective redshift, but in cosmological analyses it is also possible to take into account higher order terms e.g.  \nˆ 𝛼 ∥ , ⊥ = 𝛼 ∥ , ⊥ ( 𝑧 eff ) + 1 2 𝑑 2 𝛼 ∥ , ⊥ 𝑑𝑧 2 ! 𝑧 eff 𝜎 2 𝑧 + ... (7)  \nwhere 𝜎 2 𝑧 is the mean square width of the redshift bin. We discuss this further, and define the effective redshift, in Appendix A. Finally, in our modeling of the BAO signal we have ignored various largescale effects including wide-angle effects, unequal time effects and relativistic effects. These typically scale as ( 𝐻 / 𝑘 ) 2 in the power spectrum and so are extremely suppressed on the range of scales we fit, since our fiducial setup has 𝑘 min = 0 . 02 ℎ Mpc -1 . Many of these effects can produce out-of-phase contributions but their effect on the inferred BAO scale is rather limited (typically less than 0 . 05%; Appendix B). As such, we will not further consider them in the rest of this work.', '3 NONLINEAR BARYON ACOUSTIC OSCILLATIONS IN GALAXY CLUSTERING': 'Nonlinear structure formation and galaxy bias significantly modify the observed shape and phase of the linear BAO signal, especially pre-reconstruction (Bharadwaj 1996; Meiksin et al. 1999; Eisenstein et al. 2007a; Seo et al. 2008; Smith et al. 2008; Crocce & Scoccimarro 2008; Matsubara 2008; Padmanabhan & White 2009; Sherwin &Zaldarriaga 2012; Carlson et al. 2013). Broadly speaking, we can split the effects of these nonlinearities into two categories: (1) the damping, or modulation, of the amplitude of the linear wiggles and (2) the shift in the observed BAO scale through out-of-phase contributions to galaxy clustering.  \nThe most important contributions to both of these effects come from the displacements Ψ of galaxies on large scales, which are welldescribed by perturbation theory. Within Lagrangian perturbation theory (LPT), which models structure formation via the trajectories of galaxies, 𝒙 ( 𝒒 , 𝑡 ) = 𝒒 + Ψ ( 𝒒 , 𝑡 ) where 𝒒 is the initial Lagrangian position of the galaxy, the galaxy power spectrum in real space can be written as (Matsubara 2008; Carlson et al. 2013)  \n𝑃 ( 𝑘 ) = ∫ 𝑑 3 𝒒 𝑒 𝑖 k · 𝒒 -1 2 𝑘 𝑖 𝑘 𝑗 𝐴 𝑖 𝑗 ( 𝒒 ) GLYPH<0> 1 + 2 𝑏 𝐿 1 𝑖𝑘 𝑖 ⟨ 𝛿 ( 𝒒 ) Δ 𝑖 ⟩ + ( 𝑏 𝐿 1 ) 2 𝜉 lin ( 𝒒 ) + 2 𝑏 𝐿 1 𝑏 𝐿 2 𝑖𝑘 𝑖 ⟨ 𝛿 ( 𝒒 ) Δ 𝑖 ⟩ 𝜉 lin ( 𝒒 ) + ... GLYPH<1> (8)  \nwhere Δ = Ψ ( 𝒒 ) -Ψ ( 0 ) is the pairwise displacement between two galaxies separated by 𝒒 , 𝐴 𝑖 𝑗 = ⟨ Δ 𝑖 Δ 𝑗 ⟩ is its variance, and 𝜉 lin is the linear matter correlation function. The numbers 𝑏 𝐿 𝑛 are the free parameters in the Lagrangian bias expansion of galaxies, related to the canonical Eulerian ones through simple linear transformations e.g. 𝑏 1 = 1 + 𝑏 𝐿 1 (Desjacques et al. 2018). Throughout this work we  \nwill operate under the assumption that galaxies are biased tracers of cold dark matter and baryons but not (massive) neutrinos (Castorina et al. 2015). This has been shown to be an excellent approximation in simulations-close to the 0 . 1% level for the galaxy power spectrum in redshift space on the scales we are interested in-and as such we will operate within this regime throughout the text (Bayer et al. 2022). We note that the discussion in this section is not meant to constitute an exhaustive treatment of galaxy clustering in perturbation theory, and nor is Equation 8 meant to be the complete expression for the 1loop galaxy power spectrum; the intent is rather to pick out the pieces of physics most relevant to the BAO signal in order to estimate their systematic effect on BAO measurements. 4', '3.1 Real Space': "Let us begin with modifications to the linear BAO signal in real space due to only one power of the baryon fraction, i.e. which involves contributions of order O( 𝑃 𝑤 × 𝑃 𝑛 𝑛𝑤 ) . Physically we can think of these as distortions to the BAO feature due to gravitational nonlinearities of the smooth component.  \nWe can see the origin of both the 'in-phase' BAO damping and 'out-of-phase BAO' shift directly in Equation 8. In order to do so it is again useful to split the linear correlation function 𝜉 lin into its contributions from the wiggle and no-wiggle components 𝜉 𝑤 ( 𝑛𝑤 ) -the former has a sharp peak at 𝑟 𝑑 . Focusing on the 𝑏 2 1 piece we can see that the wiggle-component can therefore be approximated as (Vlah et al. 2016, see also e.g. Baldauf et al. 2015; Senatore & Zaldarriaga 2015; Blas et al. 2016 for discussions outside of LPT)  \n( 𝑏 𝐿 1 ) 2 ∫ 𝑑 3 𝒒 𝑒 𝑖 k · 𝒒 -1 2 𝑘 𝑖 𝑘 𝑗 𝐴 𝑖 𝑗 ( 𝒒 ) 𝜉 𝑤 ( 𝒒 ) ≈ ( 𝑏 𝐿 1 ) 2 𝑒 -1 2 𝑘 2 Σ 2 NL ∫ 𝑑 3 𝒒 𝑒 𝑖 k · 𝒒 𝜉 𝑤 ( 𝒒 ) = 𝑏 2 1 𝑒 -1 2 𝑘 2 Σ 2 NL 𝑃 𝑤 ( 𝑘 )  \nwhere Σ 2 NL is the angular average of 𝐴 𝑖 𝑗 ( 𝒒 ) for 𝑞 = 𝑟 𝑑 . At first-order in LPT, i.e. in the Zeldovich approximation, this is given by 5  \nΣ 2 NL = 2 3 ∫ 𝑑𝑘 2 𝜋 2 𝑃 ( 𝑘 ) [ 1 -𝑗 0 ( 𝑘𝑟 𝑑 )] . (9)  \nThis gives us the damping of the BAO by the bulk displacements of galaxies. Similar reasoning shows that any 'wiggle' component of the power spectrum is similarly damped (Beutler et al. 2019a; Vasudevan et al. 2019; Chen et al. 2020). An important feature of Equation 9 is that the damping comes from modes with wavelengths smaller that the BAO, i.e. with 𝑘 ≳ 1 / 𝑟 𝑑 ; physically, this reflects that galaxies separated by 𝑟 𝑑 are moved coherently by modes larger than the BAO scale, but incoherently by those smaller. In Λ CDM universes similar to our own the square displacement Σ 2 NL is dominated by linear displacements close to the 1 / 𝑟 𝑑 cutoff, whose effects  \n- 4 See Appendix C and onwards for a more complete accounting for galaxies in redshift space and post-reconstruction.  \n5 The saddle-point approximation here involves an angular integral dependent on the relative orientations of ˆ 𝑞 and ˆ 𝑘 . The approximation we have made here is to only resum the isotropic component, but other resummation choices also exist in the literature, e.g. Baldauf et al. (2015); Blas et al. (2016) have an integrand proportional to 1 -𝑗 0 + 2 𝑗 2 instead, corresponding to resumming the ˆ 𝑞 = ˆ 𝑘 component. The differences in the BAO amplitude in 1-loop perturbation theory using these two resummation choices are well below the percent level (see e.g. Figs. 3 and 5 of Chen et al. 2020) while in the BAO fitting only case their difference should be well within the priors of Σ 2 NL set in the analysis (see Section 7.2). For a more detailed exploration of the angular integral in LPT see Chen et al. (2024).  \nare sufficiently large that they have to be kept in the exponential, or resummed, instead of being perturbatively expanded. However, we stress that while BAO damping in Λ CDM universes similar to our own is dominated by contributions from relatively long-wavelength modes as described above, small-scale, nonlinear, displacements can also contribute. We discuss the importance of these effects when fitting BAO, especially considering redshift-space distortions and reconstruction, further below.  \nThe BAO shift effect can be similarly derived by approximating the nonlinear effect of no-wiggle contributions by their values near 𝑞 = 𝑟 𝑑 when they multiply well-localized wiggly components such as 𝜉 𝑤 . Let us consider the 𝑏 1 𝑏 2 term in Equation 8. In this case, writing to linear order  \n⟨ 𝛿 ( 𝒒 ) Δ 𝑖 ⟩ = -1 3 𝑞 𝑖 𝜎 2 ( 𝑞 ) (10)  \nwhere 𝜎 2 ( 𝑞 ) is the mean-square linear density fluctuation in a sphere of radius 𝑞 , 6 we have  \n2 𝑖𝑏 𝐿 1 𝑏 𝐿 2 𝑖𝑘 𝑖 ∫ 𝑑 3 𝒒 𝑒 𝑖 k · 𝒒 ⟨ 𝛿 ( 𝒒 ) Δ 𝑖 ⟩ 𝜉 𝑤 ( 𝒒 ) ≈ -2 3 𝑏 𝐿 1 𝑏 𝐿 2 𝜎 2 𝑠 𝑘 𝑖 ∇ 𝑘 𝑖 ∫ 𝑑 3 𝒒 𝑒 𝑖 k · 𝒒 𝜉 𝑤 ( 𝒒 ) = -2 3 𝜎 2 𝑠 𝑏 𝐿 1 𝑏 𝐿 2 𝑑𝑃 𝑤 𝑑 ln 𝑘 . (11)  \nIn the above we have defined 𝜎 2 𝑠 = 𝜎 2 ( 𝑟 𝑑 ) to be the mean square density at the peak of the wiggly correlation function. This directly translates into a physical shift in 𝛼 iso and, sign-included, has the interpretation that galaxies with positive 𝑏 2 are more likely to form in overdense regions, and galaxies sitting in an overdensity will tend to flowtowardseachother,shrinking the average observed BAO (Crocce &Scoccimarro 2008; Padmanabhan & White 2009; Sherwin & Zaldarriaga 2012). Importantly, this effect is caused by overdensities due to modes with wavelengths longer than 𝑟 𝑑 . While we have focused on the phase shift due to 𝑏 2 , all other quadratic contributions to clustering will produce similar shifts; working in Eulerian perturbation theory (EPT), Sherwin & Zaldarriaga (2012) showed that the real-space shift for a biased tracer with quadratic density bias is 7  \nΔ 𝛼 iso , NL = 47 105 GLYPH<18> 1 + 70 47 𝑏 2 𝑏 1 GLYPH<19> 𝜎 2 𝑠 ≈ 0 . 3% GLYPH<18> 1 + 70 47 𝑏 2 𝑏 1 GLYPH<19> GLYPH<18> 𝐷 ( 𝑧 ) 𝐷 ( 0 ) GLYPH<19> 2 (12)  \n6 This correlator, ⟨ 𝛿 Δ 𝑖 ⟩ , describes the mean inflow of matter towards an overdensity and is more commonly written in the LPT literature as (Carlson et al. 2013)  \n𝑈 𝑖 ( 𝒒 ) = -ˆ 𝑞 𝑖 ∫ 𝑑𝑘 𝑘 2 𝜋 2 𝑃 lin ( 𝑘 ) 𝑗 1 ( 𝑘𝑞 ) .  \nTo relate the two expressions we note that  \n𝑗 1 ( 𝑘𝑞 ) = 1 3 ( 𝑘𝑞 ) 𝑊 𝑞 ( 𝑘 )  \nwhere 𝑊 𝑞 = 3 ( sin ( 𝑘𝑞 )/( 𝑘𝑞 ) -cos ( 𝑘𝑞 ) )/( 𝑘𝑞 ) 2 is the Fourier transform of the spherical top hat window. This is a consequence of Birkhoff's theorem, i.e. that the inflow velocity depends only on the enclosed mass at radius 𝑞 .  \n7 Wenote that Sherwin & Zaldarriaga (2012) presented two different expressions for the BAO phase shift, the one we have quoted here and a larger one in the correlation function. However, the two expressions are in fact equivalent under Fourier transform and, since the covariance is (more) diagonal in Fourier space and therefore corrections therein more directly translated into systematic shifts in parameters, we use the former expression throughout this text. In addition, we note that the precise weighting of the 𝑘 modes in 𝜎 2 𝑠 stems from our LPT calculation.  \nFigure 1. A summary of oscillatory, non-linear, effects that lead to shifts in the measured BAO position in the real-space 1-loop matter power spectrum. Exact expressions are given in Eqs. 11-13. For reference, the dashed line is the linear theory wiggle power spectrum 𝑃 𝑤 scaled down to 2% of its nominal amplitude, and the 𝑓 2 𝑏 contribution amplified by a factor of 10, for ease of comparison.  \n<!-- image -->  \nwhere 𝐷 ( 𝑧 ) is the linear growth factor and we have evaluated the shift in the fiducial DESI cosmology.  \nThe BAO feature can also couple with itself to produce additional out-of-phase contributions to the signal. These are contributions of the order O( 𝑃 𝑛 𝑤 ) and hence 𝑓 𝑛 𝑏 . For example, at quadratic order we can write the contribution in Fourier space  \n𝑑 𝑑𝑓 2 𝑏 𝑃 ( 𝑘 ) = ∫ 𝑑 3 p ( 2 𝜋 ) 3 𝐾 2 ( p , k -p ) 𝑃 𝑤 ( 𝑝 ) 𝑃 𝑤 (| k -p |) (13)  \nwhere 𝐾 2 is the quadratic perturbation theory kernel for galaxy clustering. Since each power of 𝑃 𝑤 is sinusoidal with frequency 𝑟 𝑑 , their convolution in Fourier space can produce a frequency dependence distinct from that in linear theory. In the case of nonlinear matter clustering or linear bias we have that 𝐾 2 is the quadratic matter kernel 𝐹 2 ( p 1 , p 2 ) which peaks when p 1 = p 2 = k / 2, producing an outof-phase oscillatory signal sin ( 𝑟 𝑑 𝑘 / 2 ) 2 ∼ cos ( 𝑟 𝑑 𝑘 ) (Padmanabhan & White 2009). Since this signal only has support within one BAO wiggle of k / 2 it scales roughly as 𝑓 2 𝑏 ( 𝜋 / 𝑟 𝑑 ) 2 𝑘𝑃 𝑤 ( 𝑘 / 2 ) 2 , compared to the nonlinear shift derived from Eqs. 11 and 12 which scales as ∝ 𝜎 2 𝑠 × 𝑘𝑃 ' 𝑤 -itis thus subdominant since 𝑘 / 2 will typically appear after the peak of the power spectrum and because it contains an extra power of the baryon fraction 𝑓 𝑏 .  \nA comparison of both of the BAO shift contributions we have discussed, in the case of matter clustering, is shown in Figure 1. Note that unlike in the case of a pure sine wave the actual BAO wiggles in the linear power spectrum lead to a 𝑑 2 𝑃 / 𝑑𝑓 2 𝑏 that also contains a significant in-phase component, further suppressing the contribution of this effect to the total BAO shift.  \nIn addition to the above effects, which describe the clustering of galaxies assuming dark matter and baryons can be treated as a single matter fluid on cosmological scales, the very same physics that generates the BAO produces residual differences in dark matter and baryon clustering that remain even after decoupling. The dominant effects can be written as (Schmidt 2016; Chen et al. 2019a)  \n𝛿 𝑔 ∋ 𝑏 𝛿 𝑏𝑐 𝛿 𝑏𝑐 + 𝑏 𝜃 𝑏𝑐 𝜃 𝑏𝑐 + 𝑏 𝑣 2 𝑣 2 𝑏𝑐 (14)  \nwhere we have defined the relative linear perturbations as 𝑋 𝑏𝑐 =  \n𝑋 𝑏 -𝑋 𝑐 , with 𝛿 and 𝜃 the linear densities and velocity divergences of the baryons and dark matter, respectively. The linear relative densities and biases are constant and decay with redshift as 𝑎 -1 ( 𝑧 ) , respectively. The third term is the so-called 'relative velocity effect' which was first described in Tseliakhovich & Hirata (2010) as the effect of supersonic baryon-dark matter velocities on gravitational collapse. Unlike the first two terms, it appears only at one-loop order in perturbation theory, but can be enhanced as the relative velocity squared is a dimensionful parameter - studies of the relative velocity effect have often used the value 𝑏 𝑣 2 ∼ 0 . 01 𝜎 -2 𝑏𝑐 ( 𝑧 ) where 𝜎 2 𝑏𝑐 is the mean square fluctation of relative velocity perturbations, which correspond to roughly 30 km s -1 flows at recombination, based on early investigations into this effect.  \nCompared to the usual galaxy biases, 𝑏 1 , 𝑏 2 etc., the sizes of these relative bias parameters are somewhat less well understood, though some recent progress has been made. Barreira et al. (2020) measure the relative density bias 𝑏 𝛿 𝑏𝑐 using the separate-universe technique in hydrodynamical simulations and find distinct trends for dark matter halos (with halo mass) and galaxies (with stellar mass). While the trend in halos can be well-modelled by the change in amplitude of local fluctuations due to baryon density, as can halo-mass selected galaxies, the trend in galaxies selected by stellar mass is sensitive to other baryonic physics and exhibits an opposite (positive) trend with mass, particularly at the redshifts 𝑧 < 2 that DESI will observe. Nonetheless for all of these competing effects, the inferred bias values for both relevant halos and galaxies are well described by | 𝑏 𝛿 𝑏𝑐 | ≲ 1 for the galaxies relevant to DESI. The linear relative velocity divergence bias is less well understood but physical estimates through excursion sets and other methods (Schmidt 2016; Blazek et al. 2016) suggest | 𝑏 𝜃 𝑏𝑐 | ≲ 7 ( 1 + 𝑧 ) -1 , where the redshift dependence cancels the decay of peculiar velocities. The size of the relative velocity effect at DESI redshifts is fairly poorly understood, with estimates as high as 𝑏 𝑣 2 ∼ 0 . 01 𝜎 -2 𝑏𝑐 or as low as 10 -5 𝜎 -2 𝑏𝑐 for DESI-like galaxies (Yoo et al. 2011; Blazek et al. 2016; Schmidt 2016). A direct estimate of the effect of streaming velocities on baryonic infall into halos gives an order ( 10 5 𝑀 ⊙ / 𝑀 ) 2 / 3 effect (Blazek et al. 2016), which yields 𝑏 𝑣 2 much closer to the lower estimate for DESI galaxies, while reaching the higher value often used in the literature requires that galaxies today retain significant memory of star formation at early times (Yoo &Seljak 2013).  \nFigure 2 shows the leading contributions to the galaxy power spectrum of each of these terms in real space at 𝑧 = 0 (Yoo et al. 2011; Schmidt 2016; Blazek et al. 2016; Chen et al. 2019a), compared to linear theory wiggles. When the relative velocity effect ( 𝑏 𝑣 2 ; blue, orange and green curves) is the size of the most generous estimates quoted above it is the leading oscillatory contribution, with anti-correlated wiggles up to one quarter the size of the matter oscillations . The leading relative density contribution ( 𝑏 𝛿 𝑏𝑐 ; red) is also quite substantial, and its importance grows as a function of redshift since 𝛿 𝑏𝑐 is stationary while 𝛿 𝑚 decays as 𝐷 ( 𝑧 ) with redshift. The relative velocity divergence ( 𝑏 𝜃 𝑏𝑐 ; purple) on the other hand is suppressed by up to two orders of magnitude compared to these other contributions and thus its contribution to any BAO shift should be negligible. It is worth noting that all of the shown relative bias contributions contain significant or even dominant in-, or exactly out-of-phase contributions with oscillations with the opposite sign to the BAO. These oscillations mainly modify the amplitude ofthe BAO rather than shift it; if these are present at a significant level in DESI galaxies they would lead to modifications of their clustering (e.g. inconsistencies between the BAO and broadband amplitudes) which could be constrained in full-shape fits (Blazek et al. 2016;  \nFigure 2. Contributions to the real-space galaxy power spectrum at 𝑧 = 0 from relative baryon-dark matter perturbations, compared to the matter 𝑃 𝑤 (black). All three leading effects show prominent oscillatory features independent from the 𝑃 𝑤 , with the largest contribution from the relative density perturbation for values 𝑏 𝛿 𝑏𝑐 ∼ 1 (red). The contributions due to the relative velocity divergence ( 𝑏 𝜃 𝑏𝑐 ) (purple) are quite small and are scaled here by a factor of 100 for visibility. The relative velocity effect terms ( ∝ 𝑏 𝑣 2 ) assume an effect of order 0 . 01 𝜎 -2 𝑣 , a very pessimistic choice at the upper end of the range of estimates. Contributions due to the cross correlation of the 𝑏 𝑣 2 term with matter, linear and quadratic bias are shown in blue, orange and green, respectively, with the former two showing particularly prominent oscillations at this value of 𝑏 𝑣 2 .  \n<!-- image -->  \nBeutler et al. 2017b; Slepian et al. 2018). These constraints could then be used to reduce the theoretical error on the BAO.  \nOverall, the unknown amplitude of these effects means they must be accounted for, conservatively, in our theoretical systematic error budget. However, this requires further consideration of these effects in redshift-space, which will be carried out next.", '3.2 Redshift Space': "We now extend the discussion in the previous section to include redshift-space distortions. In LPT, redshift-space distortions can be modelled by including an extra displacement along the line of sight ˆ 𝑛 proportional to the peculiar velocity, i.e. Ψ 𝑠 = Ψ + ˆ 𝑛 GLYPH<0> ˆ 𝑛 · / Ψ GLYPH<1> /( 𝑎𝐻 ) (Matsubara 2008; Carlson et al. 2013). Within linear theory this is equivalent to taking Ψ 𝑠,𝑖 = 𝑅 𝑖 𝑗 Ψ 𝑗 , where 𝑅 𝑖 𝑗 = 𝛿 𝑖 𝑗 + 𝑓 ˆ 𝑛 𝑖 ˆ 𝑛 𝑗 and 𝑓 ( 𝑧 ) is the linear growth rate.  \nThe effect of redshift-space distortions on the BAO damping is straightforward. Since the exponentiated displacements all get multiplied by 𝑅 𝑖 𝑗 , the damping factor acquires an anisotropy given by (Eisenstein et al. 2007a; Matsubara 2008; Ivanov & Sibiryakov 2018)  \n-1 2 𝑘 𝑖 𝑘 𝑗 𝑅 𝑖𝑎 𝑅 𝑗𝑏 Σ 2 NL 𝛿 𝑎𝑏 = -1 2 𝐾 2 Σ 2 NL , 𝐾 𝑖 = 𝑅 𝑖 𝑗 𝑘 𝑗 , (15)  \nsuch that 𝐾 2 = ( 1 + 𝑓 ( 2 + 𝑓 ) 𝜇 2 ) 𝑘 2 . As in real space, small-scale nonlinear effects can also affect the damping of the BAO peak. In practice, the most significant contribution will be from the so-called 'Fingers of God' (FoG) - random motions, and therefore redshiftspace displacements - due to virial motions of galaxies within halos (Jackson 1972). To leading order, these contribute an additional -1 2 𝑘 2 𝜇 2 Σ 2 FoG to the exponential damping. Since virial velocities  \ncome from a range of halo masses with different dispersions (Sheth 1996; Diaferio & Geller 1996; White 2001) their distribution is not Gaussian (Peebles 1976; Davis & Peebles 1983). Their overall effect may be better described by other probability distributions (Seljak & McDonald 2011; Chen et al. 2021), such as the e.g. Lorentzian form sometimes used in the literature (Peacock 1992; Park et al. 1994; Peacock & Dodds 1994; Ballinger et al. 1996). We will quantify the impact of these choices in §7.2.  \nThe effect on the nonlinear BAO shift is somewhat more complicated. Following the same reasoning as in real space, we can easily come to the conclusion that the BAO shift should become anisotropic and enhanced along the line of sight. This is because redshift-space distortions from peculiar velocities streaming towards overdensities will cause galaxies separated by an overdensity along the line of sight to appear closer to each other than those separated perpendicular to such an overdensity. An exact calculation of all the anisotropic nonlinearities proportional to the large-scale overdensity 𝜎 2 𝑠 yields (Appendix C)  \n𝑃 shift = -( 1 / 315 ) 𝜎 2 𝑠 𝑘𝑃 ' 𝑤 ( 𝑘 )( 𝑏 1 + 𝑓 𝜇 2 )( 141 𝑏 1 + 210 𝑏 2 + 112 𝑏 𝑠 + 𝑓 ( 93 + 330 𝑏 1 + 210 𝑏 2 + 112 𝑏 𝑠 -42 ( 1 + 𝑏 1 ) 𝑓 ) 𝜇 2 + 3 ( 80 + 77 𝑏 1 -28 𝑓 ) 𝑓 2 𝜇 4 + 189 𝑓 3 𝜇 6 ) , (16)  \nwhere we have also included the possible effect of a tidal bias 𝑏 𝑠 . Evidently, the shift can no longer be mapped into a single Δ 𝛼 iso as it could in real space, though it is worth noting that the total shift for modes perpendicular to the line of sight ( 𝜇 = 0) remains the same, while the shift along the line of sight becomes significantly enhanced.  \nIn order to translate Equation 16 into systematic shifts in 𝛼 ∥ , 𝛼 ⊥ we conduct a simple Fisher analysis assuming a cosmic-variance limited survey and Gaussian covariances with a fiducial model given by the Kaiser form plus a polynomial broadband model as described in §5. Defining the systematic error vector 𝜖 = ( 𝑃 0 shift , 𝑃 2 shift ) to be the multipoles of Equation 16 we have that the shift in 𝛼 's is given by  \nΔ 𝛼 NL ∥ , ⊥ = 𝐹 -1 𝛼 ∥ , ⊥ ,𝑏 𝑡 𝑏 𝑖 𝐶 -1 𝑖 𝑗 𝜖 𝑗 , (17)  \nwhere 𝐶 -1 𝑖 𝑗 is the inverse covariance, 𝑡 𝑏 𝑖 are the linear templates for the free parameters 𝑝 𝑏 and 𝐹 𝑎𝑏 is the Fisher matrix. For simplicity we assume the coevolution relations 𝑏 2 = 𝑏 𝐿 2 + 8 21 𝑏 𝐿 1 and 𝑏 𝑠 = -2 7 𝑏 𝐿 1 , where the Lagrangian 𝑏 𝐿 2 is fixed to be a function of 𝑏 𝐿 1 according to the Sheth-Tormen mass function (Sheth et al. 2001).  \nFigure 3 shows the expected shifts in the BAO scale in the above setup as a function of 𝑏 1 and redshift. The solid lines show the expected BAO shift in the absence of any nonlinear bias: they roughly reproduce the real-space prediction (black dotted) in the case of 𝛼 ⊥ but are up to 3 to 4 times larger for 𝛼 ∥ . This makes sense since the perpendicular part of Equation 16 is identical to the real-space prediction.  \nThe dashed lines in Figure 3 show the result of adding in nonzero quadratic density and shear biases using peak-background split predictions. These show significant deviations away from the linear bias-only predictions; roughly, for 𝑏 𝐿 1 < 1the quadratic bias becomes negative, reducing the effect, and vice versa for more biased tracers, which cluster more strongly around overdensities. The bias relations used here are known to qualitatively reproduce the behaviour of dark matter halos in N-body simulations (Lazeyras & Schmidt 2018; Abidi & Baldauf 2018); however, they are not guaranteed to hold for the galaxy samples observed in DESI, and as such the scatter between the different lines in the plot should be thought of as an indicator of the  \nFigure 3. Fisher forecasts for the BAO shift in redshift space due to nonlinear clustering as a function of redshift 𝑧 and linear bias 𝑏 1 , assuming a signaldominated survey. Solid lines show results for linear bias only while dashed lines show the shift when a quadratic bias parameter, predicted from the peakbackground split, is included, showing that the scatter in potential shifts due to the unknown nonlinear bias parameters is substantial. The dotted black line shows the shift in real space, which is similar in size to Δ 𝛼 ⊥ in most cases but substantially smaller than Δ 𝛼 ∥ .  \n<!-- image -->  \npossible systematic scatter due to the unknown size of nonlinear biasing. Conversely, a full-shape fit of galaxy clustering aimed not just at extracting the BAO signal and where all perturbation theory parameters are varied automatically incorporates this uncertainty into its parameter constraints. At least in the case of pre-reconstruction, where perturbation theory models are quite mature, an alternative procedure to fitting the BAO may be to employ such a model, as one would do in a 'full' fit of the galaxy power spectrum, while simultaneously marginalizing over broadband terms (§2) as is typically done in BAO fits. Beyond standard BAO measurements, we note that these shifts have a larger effect on 𝛼 ∥ , ⊥ than the expected shift due to light species such as neutrinos, which has emerged as a topic of recent interest (e.g. Baumann et al. (2019)), and therefore must be carefully marginalized over when looking for 'exotic' effects in BAO physics; on the other hand, as we will see, they are very well mitigated by reconstruction.  \nConverting the error estimates above into 𝛼 iso and 𝛼 ap , the former unsurprisingly is subject to nonlinear shifts spanning the geometric average Δ 𝛼 ∥ / 3 + 2 Δ 𝛼 ⊥ / 3, such that the range of scatter due to the unknown size of nonlinear biases is directly translated over from Figure 3. In the case that the nonlinear biases are given by the peak background split and Lagrangian coevolution we have, roughly, Δ 𝛼 iso ≈ 0 . 5 𝜎 2 𝑠 ( 1 + ( 𝑏 1 -1 )) 𝐷 ( 𝑧 ) 2 , or roughly 0 . 3% for a bias one tracer at redshift zero, while neglecting any nonlinear  \nEulerian bias gives ≈ 0 . 9 𝜎 2 𝑠 𝐷 ( 𝑧 ) 2 . The scatter between these two again is indicative of the uncertainty induced in pre-reconstruction measurements by nonlinearities. On the other hand, it is interesting to note that this scatter essentially cancels in 𝛼 ap such that we have Δ 𝛼 ap ≈ 1 . 2 𝜎 2 𝑠 𝑓 ( 𝑧 ) 𝐷 ( 𝑧 ) 2 , corresponding to approximately a 0 . 4% shift at 𝑧 = 0, to within 0 . 1% accuracy over the entire range of biases and redshifts shown, with only minor variations increasing with nonlinear bias. This suggests that a direct correction based on this estimate can be performed for 𝛼 ap when measuring it from pre-reconstruction BAO. We caution that the fitting forms above are approximate and intended for error estimation.  \nWe can also apply the same Fisher formalism to extract the BAO shifts due to the quadratic coupling of 𝑃 𝑤 with itself and the relative velocity effect. Redshift-space distortions similarly enhance the size of the former along the line of sight; however, since this effect ultimately has a rather distinguishable scale dependence compared to the BAO wiggles and their derivatives the contribution to the measured BAO scale is extremely suppressed - for a bias 𝑏 = 2 tracer at 𝑧 = 0 . 8 the effect is of order 𝑓 2 𝑏 × 0 . 1%, which is entirely negligible. This result justifies neglecting any order 𝑓 2 𝑏 effects in BAO fits even though they are nominally only suppressed by one power of the baryon fraction 𝑓 𝑏 ≈ 15%.  \nFinally, we compute the redshift-space power spectrum contributions of the relative velocity effect using the Lagrangian formalism of Chen et al. (2019a), converting their real-space expressions into redshift-space as in Vlah & White (2019); Chen et al. (2021) (see also Schmidt 2016). By operating within LPT we avoid having to perform wiggle/no-wiggle splits on the relative velocity spectra to model their nonlinear damping; the standard algorithms to perform this split are optimized for the matter power spectrum and cannot typically be used 'out of the box' in this case where the wiggles are also more pronounced. Keeping corrections only to first order in the relative bias parameters we can express the resulting BAO shifts from Equation 17 as Δ 𝛼 ∥ , ⊥ = 𝑏 𝑖 ( 𝑑𝛼 ∥ , ⊥ / 𝑑𝑏 𝑖 ) for 𝑏 𝑖 = { 𝑏 𝛿 𝑏𝑐 , 𝑏 𝜃 𝑏𝑐 , 𝑏 𝑣 2 } . In addition to the shift in best-fit parameters, the presence of these relative bias terms can also change the strength of BAO constraints by modifying the BAO signal amplitude with their in-phase components. However, this effect will be avoided so long as the covariance matrices and model can match the amplitude of the observed BAO.  \nFigure 4 shows the resulting redshift-space BAO shift per unit relative bias. The top panel shows the effect of the relative overdensity bias 𝛿 𝑏𝑐 . The magnitude of the shift diminishes towards lower redshift and higher biases and is opposite in sign for 𝛼 iso , ap . This is as expected, since the relative component stays constant with redshift unlike the growing adiabatic component, so perpendicular to the line of sight the relative size of the effect scales as 𝑏 𝛿 𝑏𝑐 /( 𝑏 1 𝐷 ( 𝑧 )) , while along the line of sight, the shift is somewhat reduced by the presence of redshift-space distortions sourced by matter, such that the shift in 𝛼 ap is negative. We note that while Barreira et al. (2020) find that 𝑏 𝛿 𝑏𝑐 is typically of order one for both galaxies and halos and the expected shifts are therefore of the same order of magnitude as the derivatives shown in the top panel, they also find that 𝑏 𝛿 𝑏𝑐 tends to increase with halo mass. Given that the exact value of this bias is unknown as a conservative estimate we should therefore expect a ∼ 0 . 05 𝐷 -1 ( 𝑧 ) %shift in the BAO from this effect.  \nThe bottom panel similarly shows the expected shift per 0 . 01 𝑏 𝑣 2 as a function of bias and redshift. Since these terms exist at 1-loop order it is necessary to make some assumptions about nonlinear biasing; wehave therefore adopted the peak-background split values of 𝑏 2 like in Figure 3, though our results in this case are not very sensitive to this choice as the relative velocity effect is mostly sourced by long-  \nFigure 4. BAO shifts from relative velocity effect as a function of linear bias 𝑏 1 and redshift 𝑧 , computed for (top) the relative density bias 𝑏 𝛿 𝑏𝑐 = 1 and (bottom) 𝑏 𝑣 2 = 0 . 01 𝜎 -2 𝑏𝑐 . The normalization for 𝑏 𝑣 2 here is on the upper range of estimates in the literature, with estimates in the literature spanning three orders of magnitude smaller-we adopt 10 -3 𝜎 -2 𝑏𝑐 as the fiducial mean in this work. We have omitted the shift due to the relative velocity divergence since its effect is expected to be much smaller than the two shown.  \n<!-- image -->  \nh modes. Here again the signal decreases as a function of linear bias, but unlike the relative density bias, it stays roughly constant with redshift since the 1-loop terms contributing to it scale similarly to the linear adiabatic modes. Given the lack of a clear measurement in either data or simulations of the value of 𝑏 𝑣 2 , or a clear theoretical consensus in the literature about its expected size, the derivatives plotted can be thought of as a worst-case scenario for this bias. If we take a log-uniform prior on 𝑏 𝑣 2 spanning the range 10 -5 to 10 -2 𝜎 -2 𝑏𝑐 the mean shift is roughly 0 . 02 𝑏 -1 1 % on 𝛼 iso with negligible redshift dependence and a significantly smaller anisotropic shift. However, since the straightforward estimate based on halo mass yields of a few times 10 -5 , absent discovering an astrophysical mechanism that preserves the relative velocity signal in early star formation in DESI target galaxies it seems safe to assume this is also a conservative estimate. Further, while the shift due to a 𝑏 𝑣 2 in the upper range of reported values in the literature would be comparable to the statistical precision of the full Y5 DESI data, we expect that full-shape analyses in the meantime will be able to place stronger bounds on its value, since data from the BOSS survey already come quite close to constraining these values (Yoo & Seljak 2013; Beutler et al. 2017b).", '4.1 The Standard Reconstruction Algorithm': "In §3, we considered the effects of unadulterated nonlinearities in structure formation on the BAO signal in galaxy clustering. However, as described in the introduction (§1), the bread-and-butter of BAOmeasurements in spectroscopic surveys like DESI concerns the BAOsignal in galaxy clustering post-reconstruction . The reconstruction algorithm presented in Eisenstein et al. (2007b) amounts to a nonlinear transformation of the observed galaxy density field to enhance and unbias the BAO signal - the purpose of this subsection is to review these effects in some detail using perturbation theory (see Padmanabhan et al. 2009; Noh et al. 2009; White 2015; Schmittfull et al. 2015; Hikage et al. 2017; Ding et al. 2018; Chen et al. 2019b; Hikage et al. 2020, for earlier discussions).  \nThe first step of the reconstruction algorithm is to generate the reconstructed displacement Ψ 𝑟 ( s ) at a given redshift space position s . Schematically, this is done by dividing the galaxy density field by 𝑏 + 𝑓 𝜇 2 after smoothing the observed galaxy density field by a Gaussian filter S( 𝑘 ) = 𝑒 -𝑘 2 𝑅 2 / 2 with 𝑅 sufficiently large that the smoothed field is well-described by linear theory. In this case, the linear Zeldovich displacement can be reconstructed as  \nΨ 𝑟 ( k ) = 𝑖 k 𝑘 2 GLYPH<18> S( 𝑘 ) 𝛿 𝑔 ( k ) 𝑏 + 𝑓 𝜇 2 GLYPH<19> ≈ Ψ Zel ( k ) . (18)  \nIn practice, survey geometry and other effects mean that the simple procedure above cannot be followed literally (i.e. the un-windowed 𝛿 𝑔 ( k ) cannot be directly observed and divided by the Kaiser factor). In the below, however, we will mostly use the approximation in Equation 18 except when describing potential biases due to imperfect reconstruction, leaving the details of the reconstruction algorithm itself to our companion 'optimal reconstruction' papers (Chen 2024; Paillas et al. 2024).  \nGiven Ψ 𝑟 , the reconstruction algorithm proceeds to shift both the observed galaxies and randoms by the negative reconstructed displacement in order to cancel the nonlinear damping due to longwavelength displacements, and the reconstructed density field is defined to be the difference between the displaced galaxies and shifted randoms, 𝛿 𝑟 = 𝛿 𝑑 -𝛿 𝑠 . Subtracting off the overdensity of the shifted randoms restores the linear clustering power removed by displacing the galaxies. Using number conservation, the overdensities of the galaxies and randoms post-reconstruction can be expressed as (Schmittfull et al. 2015)  \n1 + 𝛿 𝑑,𝑠 ( 𝒙 𝑟 ) = ∫ 𝑑 3 𝒙 𝛿 𝐷 ( 𝒙 𝑟 -𝒙 -Ψ 𝑟 ( 𝒙 )) GLYPH<0> 1 + 𝛿 𝑑,𝑠 ( 𝒙 ) GLYPH<1> . (19)  \nHere 𝒙 𝑟 is the position of galaxies and randoms post-reconstruction, and the densities 𝛿 𝑑,𝑠 ( 𝒙 ) refer to their pre-reconstruction densities, i.e. 𝛿 𝑑 ( 𝒙 ) = 𝛿 𝑔 ( 𝒙 ) and 𝛿 𝑠 ( 𝒙 ) = 0. Equivalently we have in Fourier space  \n𝛿 𝑑,𝑠 ( k ) = ∫ 𝑑 3 𝒙 𝑒 -𝑖 k · 𝒙 GLYPH<16> 𝑒 -𝑖 k · Ψ 𝑟 ( 𝒙 ) GLYPH<0> 1 + 𝛿 𝑑,𝑠 ( 𝒙 ) GLYPH<1> -1 GLYPH<17> . (20)  \nIn order to make a connection with Lagrangian perturbation theory we can further substitute 𝒙 = 𝒒 + Ψ ( 𝒒 ) , such that the integral over 𝒙 is replaced with 𝒒 and 𝒙 𝑟 = 𝒒 + Ψ ( 𝒒 ) + Ψ 𝑟 ( 𝒒 + Ψ ) (White 2015; Chen et al. 2019a).  \nIn order to include redshift-space distortions the reconstructed displacement is multiplied by the matrix 𝑅 𝑖 𝑗 (§3.2). Where this boost factor is applied varies across the literature; in this paper, and generally in DESI Y1, we will focus on two conventions:  \n- · RecIso : galaxies shifted by -R Ψ 𝑟 , randoms shifted by -Ψ 𝑟 𝛿 𝑑 = ( 𝑏 1 + 𝑓 𝜇 2 - ( 1 + 𝑓 𝜇 2 )S) 𝛿, 𝛿 𝑠 = -S 𝛿 , 𝛿 𝑟 = ( 𝑏 1 + ( 1 - S) 𝑓 𝜇 2 ) 𝛿 ,\n- · RecSym : galaxies shifted by -R Ψ 𝑟 , randoms shifted by -R Ψ 𝑟 . 𝛿 𝑑 = ( 𝑏 1 + 𝑓 𝜇 2 - ( 1 + 𝑓 𝜇 2 )S) 𝛿, 𝛿 𝑠 = -( 1 + 𝑓 𝜇 2 )S 𝛿 , 𝛿 𝑟 = ( 𝑏 1 + 𝑓 𝜇 2 ) 𝛿 ,  \nwhere we have also written down the linear theory predictions for each scheme. The former convention, RecIso , was chosen by the SDSS, BOSS and eBOSS collaborations (Padmanabhan et al. 2012; Alam et al. 2017, 2021); since linear redshift space distortions are sourced by the same linear velocities as in R Ψ 𝑟 , removing them from the displaced galaxies without doing the same to the randoms leads to the removal of linear redshift space distortions on large scales and 'isotropizes' the clustering. The second convention, RecSym , was introduced in White (2015) and treats randoms and galaxies 'symmetrically' such that the linear theory broadband is preserved. It is important to note that (a) the anisotropic clustering still contains BAO information, so subtracting it may result in weaker BAO constraints for a fixed shot noise and (b) the cancellation of RSD depends on the extent to which the true and fiducial cosmologies match, such that even given the true Zeldovich displacement the anisotropic term would be more like ( 𝑓 - S 𝑓 fid ) 𝜇 2 , leaving a residual that could be confused for anisotropic BAO from the Alcock-Paczynski effect. RecSym is on the other hand robust against any such mistakes to leading order since both particles and randoms are moved identically.", '4.2 Nonlinear Modeling of the Reconstructed Field': "The reconstructed field 𝛿 𝑟 carries a somewhat different set of nonlinearities than the original galaxy density 𝛿 𝑔 , as expected since the reconstruction algorithm was designed specifically to reduce the nonlinear damping of the BAO signal. The nonlinear modelling of the reconstructed BAO is rather nontrivial and, while much progress has been made in recent years including efforts to model the postreconstruction matter field (Hikage et al. 2017, 2020) and galaxies in the Zeldovich approximation (White 2015; Chen et al. 2019b), to the best of the authors' knowledge no complete model taking into account both nonlinear dynamics and bias as well as dependences on the reconstruction algorithm itself (e.g. due to the fiducial cosmology assumed) has yet been presented. Nonetheless, it is possible to understand the effect of nonlinearities post-reconstruction using the same methods as the pre-reconstruction case. In the below, we will outline general arguments. More detailed calculations, particularly in the more tractable case of matter, will be presented in Appendix D.", '4.2.1 Real Space': "Let us first examine the damping of the BAO peak postreconstruction. For the displaced galaxies we can write the following version of Equation 20:  \n𝛿 𝑑 ( k ) = ∫ 𝑑 3 𝒒 𝑒 𝑖 k · ( 𝒒 + Ψ ( 𝒒 )+ Ψ 𝑟 ( 𝒒 + Ψ ) ) GLYPH<0> 1 + 𝛿 𝑔 ( 𝒒 ) GLYPH<1> , (21)  \nwhere we have dropped delta functions in k for brevity. Compared to the pre-reconstruction case we can see that the featureless displacement that smooths the power spectrum is now given by Ψ ( 𝒒 ) + Ψ 𝑟 ( 𝒒 ) ≈ ( 1 -𝑆 ) ∗ Ψ to leading order, i.e. the damping parameter is now  \nΣ 2 rec = 2 3 ∫ 𝑑𝑘 2 𝜋 2 ( 1 - S) 2 𝑃 ( 𝑘 )( 1 -𝑗 0 ( 𝑘𝑟 𝑑 )) . (22)  \nThe case of the shifted field is similar: At first sight, it looks like the shifted field starts out un-displaced so that the only BAO damping comesfromtheshift we apply to the randoms; however, since the shift field is computed from the Eulerian density field it already comes pre-displaced, such that its ultimate damping scale is the same as in the displaced field (Appendix D2). Another way to see this is that the full reconstructed density field is given via Equation 20 as (see also Shirasaki et al. 2021; Sugiyama et al. 2021)  \n𝛿 rec ( k ) = ∫ 𝑑 3 𝒙 𝑒 -𝑖 k · 𝒙 -𝑖 k · Ψ 𝑟 ( 𝒙 ) 𝛿 𝑔 ( 𝒙 ) (23)  \nand, since 𝛿 𝑔 ( 𝒙 ) can be expressed as a sum of bias operators at the Lagrangian position 𝒒 (Mirbabayi et al. 2015) it follows that the total shift is given in Equation 22. This result that the BAO post-reconstruction is still damped by a simple single-exponential form, which we derive for the first time, is in contrast to previous derivations in the literature using perturbation theory and EFT Ding et al. (2018); Chen et al. (2019b).  \nWe thus have that in real space the damping form is equivalent to the pre-reconstruction case except with a modified damping parameter, i.e.  \n𝑃 ( 𝑘 ) = 𝑏 2 GLYPH<16> 𝑃 𝑛𝑤 ( 𝑘 ) + 𝑃 𝑤 ( 𝑘 ) 𝑒 -1 2 𝑘 2 Σ 2 rec GLYPH<17> . (24)  \nIn addition, since the shift field is reconstructed from noisy galaxy densities it will add an additional uncorrelated damping proportional to the shot noise 𝑃 shot / ¯ 𝑛 (White 2010)  \nΣ 2 shot = 2 3 𝑃 shot 𝑏 2 ¯ 𝑛 ∫ 𝑑𝑘 2 𝜋 2 S( 𝑘 ) 2 ( 1 -𝑗 0 ( 𝑘𝑟 𝑑 )) ≈ 1 6 𝜋 3 / 2 𝑃 shot 𝑏 2 ¯ 𝑛𝑅 , (25)  \nwhere we have assumed that 𝑅 ≪ 𝑟 𝑑 . Here 𝑃 shot describes the deviation of the galaxy stochasticity away from a Poisson distribution, where we would have 𝑃 shot = 1. For galaxies with inverse comoving density 3 × 10 3 ℎ -3 Mpc 3 this is equal to Σ 2 shot = 6 𝑏 -2 𝑅 -1 15 𝑃 shot ℎ -2 Mpc 2 where 𝑅 15 is the reconstruction smoothing in units of 15 ℎ -1 Mpc. Again, it is important to note that this will not in general be the only additional contribution to the postreconstruction damping, but its presence is a good gauge of the size of nonlinear corrections to the Zeldovich damping. Some illustrative values of the expected damping as a function and smoothing scale are shown in Figure 5 - while in general the damping can be somewhat reduced by decreasing 𝑅 , other nonlinearities eventually make this choice sub-optimal. In the case of bias values and redshifts representative of the DESI QSO sample, we can see that the extra damping incurred by the shot noise brings the total post-reconstruction damping above the pre-reconstruction one for essentially all smoothing scales. However, for a suitably large 𝑅 the smoothing scales are rather comparable while the nonlinear BAO shift is still cancelled post-reconstruction such that performing reconstruction may still be preferred.  \nThe effect of reconstruction on the nonlinear BAO shift follows a similar logic. Reconstruction explicitly removes, up to the smoothing scale 𝑅 , the Zeldovich displacement terms that produce the shift; since the shift depends on modes 𝑘 ≲ 1 / 𝑟 𝑑 ≪ 1 / 𝑅 the cancellation is essentially perfect if the linear displacement is properly estimated (see Appendix D1 for an explicit calculation). However, misestimating the linear bias 𝑏 1 can lead to residual displacements postreconstruction, leading to errors roughly of order ( Δ 𝑏 1 / 𝑏 1 ) Δ 𝛼 iso , NL , where Δ 𝑏 1 is the difference between the true linear bias and that assumed during the application of the reconstruction algorithm. It is important to note that the relevant quantity here is the error on the linear bias and not its product with 𝜎 8 . The latter is typically extremely well constrained by the data while, absent external priors,  \n<!-- image -->  \nFigure 5. Contributions to the post-reconstruction BAO damping due to Zeldovich displacements (blue), shot noise (orange) and nonlinearities (green) as a function of smoothing scale 𝑅 , with the total shown in the black dashed line. The black dotted line shows the full Zeldovich displacement at the BAO scale, i.e. the pre-reconstruction damping scale. The green curve is estimated by computing the mean square reconstructed displacement ⟨ Ψ 2 𝑟 ⟩ using the 1-loop power spectrum for the value of bias shown, with nonlinear Eulerian biases assumed to be zero. The top and bottom panels show representative dampings for the DESI LRG and QSO samples, respectively, with the shot noise value reported in [ ℎ -3 Mpc 3 ].  \n<!-- image -->  \nthe constraints on the former are essentially dictated by the errors on 𝜎 8 itself. Conceptually, this is because we can only cancel the shift as well as we can know the amplitude of matter fluctuations within spheres of the BAO radius 𝑟 𝑑 .  \nReconstruction also has a mixed effect on the relative biases. The contributions due to 𝑏 𝛿 𝑏𝑐 , 𝜃 𝑏𝑐 are essentially linear BAO's that are out-of-phase with the matter one, and reconstruction thus simply sharpens them in the same way, taking the damping from Σ 2 NL → Σ 2 rec . Since both signals are identically strengthened, reconstruction is not expected to reduce the theoretical error due to this effect. On the other hand, the 𝑏 𝑣 2 contributions are due to its coupling with quadratic nonlinearities in the galaxy density, which is in fact dominated by the large-scale displacements cancelled by reconstruction as 𝑣 2 𝑏𝑐 is dominated by long-wavelength modes (Blazek et al. 2016). While we thus expect reconstruction to reduce the contributions of the relative velocity effect to the BAO-and indeed calculations within the Zeldovich approximation show that its oscillatory components are reduced by roughly a factor of two postreconstruction-a complete calculation of the post-reconstruction relative velocity effect is outside the scope of this work. A conservative estimate is therefore to use the pre-reconstruction shifts as shown in Figure 4 in our error budget.  \nIn the case of RecSym the arguments in real space carry over to redshift space essentially identically; since both the reconstructed displacement and the Zeldovich displacements responsible for the nonlinear BAO shift and BAO damping are multiplied by the same matrix R , reconstruction simply cancels the shift equally in redshift space and we get 𝑘 2 Σ 2 rec → 𝐾 2 ( 𝜇 ) Σ 2 rec as in the pre-reconstruction case. This motivates the fitting form  \nC( 𝑘, 𝜇 ) = ( 𝑏 + 𝑓 𝜇 2 ) 2 exp GLYPH<20> -1 2 𝑘 2 GLYPH<16> 𝜇 2 Σ 2 ∥ + ( 1 -𝜇 2 ) Σ 2 ⊥ GLYPH<17> GLYPH<21> (26)  \nboth pre- and post-reconstruction, with the understanding that the phase shift on top of this will be larger pre-reconstruction, and that Σ 2 ∥ , ⊥ are roughly given by the Zeldovich prediction with uncertainty due to nonlinearities and imperfect knowledge of the true cosmology. In particular, as with the linear bias, the extent to which these cancellations happen depends on the fidelity with which we can constrain the linear growth rate 𝑓 , such that we should expect a remaining theoretical error of size ( Δ 𝑓 / 𝑓 ) Δ 𝛼 ⊥ , ∥ . Both of these errors are related to our certainty on the underlying matter clustering amplitude 𝜎 8 , which we can conservatively place a 5% uncertainty on given current discrepancies between the primary CMB and late-universe probes, though we note that differences up to 10% with the primary CMB exist e.g. from galaxy redshift-space clustering (Abdalla et al. 2022). We remind the reader that, while Zeldovich displacements form the dominant contribution to BAO damping, small-scale nonlinearities can also contribute and will tend to be equivalent to Equation 26 only at leading order. We investigate the effect of including such effects (e.g. FoG) into BAO fits in §7.2, finding that they have negligible impact even at DESI precision. More generally, these effects beyond leading order can be mitigated by limiting the range of scales fit, especially in the power spectrum.  \nIn the case of RecIso however some care must be taken, because only the real-space portion of the Zeldovich displacement is cancelled for the shifted randoms. In this case, it is reasonable to assume that much (though not all) of the contributions due to the anisotropic BAO shift will remain, though the 𝜇 = 0 contribution should still be cancelled. We show this explicitly for matter in Appendix D3. Similarly, while the displaced galaxies will still have their nonlinear damping reduced to 𝐾 2 Σ 2 rec the shifted fields and their cross spectra with galaxies will now be damped according to  \nΣ 2 𝑠𝑠 ( 𝜇 ) = Σ 2 rec + 𝑓 2 𝜇 2 Σ 2 NL + 2 𝑓 𝜇 2 Σ 2 𝑥 , Σ 2 𝑑𝑠 ( 𝜇 ) = ( 1 + 𝑓 𝜇 2 ) Σ 2 rec + 𝑓 ( 1 + 𝑓 ) 𝜇 2 Σ 2 𝑥 + 𝑓 2 𝜇 2 Σ 2 S (27)  \nwhere we have defined Σ 2 𝑥 = 2 3 ∫ 𝑑𝑘 2 𝜋 2 ( 1 - S) 𝑃 lin ( 𝑘 ) ( 1 -𝑗 0 ( 𝑟 𝑑 𝑘 )) as the cross-correlation between the pre- and post-reconstruction Zeldovich displacements and Σ 2 S = 1 3 ∫ 𝑑𝑘 2 𝜋 2 S 2 𝑃 lin ( 𝑘 ) is the meansquare Zeldovich displacement after smoothing by S . The distinct pairwise displacements above between the shifted galaxies and randoms in RecIso mean that the reconstructed power spectrum in this case does not have BAO which are damped by a single exponential, but rather piecewise according to Equation 27, which can be also seen from the fact that Equation 23 does not hold if the galaxies and randoms are not moved by the same amount. The above expressions reduce to the pre-reconstruction displacement when S → 0, and the post-reconstruction displacement in real space if 𝑓 → 0, but in general the BAO signal in RecIso will be more damped than it is in RecSym , and the nonlinear shift larger, because of the incomplete cancellation of displacements along the line of sight. We thus see that RecIso suffers from (a) a partially uncancelled nonlinear BAO  \nshift (b) more damping (c) less signal and (d) greater sensitivity to systematics in reconstruction than RecSym .", '5 BROADBAND MODELS IN FOURIER AND CONFIGURATION SPACE': "Thus far, we have developed the theory for the effects of nonlinear structure formation on the oscillatory BAO signal. The main contributions come from large-scale displacements which smoothly alter the envelope of the damped BAO wiggles and shift the observed BAO.These can be effectively removed through reconstruction, leading to a robust BAO signal that can be only weakly biased through second-order bias effects and (for pre-reconstruction or under RecSym convention) damped by a single exponential damping term. However, nonlinearities-which are modified, but not removed, by reconstruction- not to mention systematics in the data and differences in the fiducial and true cosmologies, also contribute to the broadband shape of the observed galaxy 2-point function. Since the express purpose of BAO fits is to isolate the geometric information in the oscillatory part of the power spectrum only, the ability to marginalize over any broadband differences is an integral part of any BAO modeling.  \nAnessential feature of such a broadband model is that it must fit any feature in the power spectrum with frequency in Fourier space lower than that of the BAO, thereby isolating the BAO information from the rest of the clustering signal. This is equivalent in configuration space to saying that it must remove any possible features localized at radii smaller than the BAO scale, as we will discuss in detail in §5.2. It is important to emphasize that this is done regardless of the source (fiducial cosmology, observational systematic or nonlinearity) of the broadband contribution and, conversely, if a systematic contributes to clustering in a way that is not sufficiently smooth it would need to be specifically modeled in order to mitigate its effect on the BAO.  \nThe broadband terms in our fiducial BAO model (Equation 5) are contained in B( 𝑘, 𝜇 ) 𝑃 𝑛𝑤 + D( 𝑘, 𝜇 ) . In principle D can be any sufficiently flexible fitting form non-degenerate with the BAO wiggles, though in practice it is useful to isolate a 'baseline' contribution expected from linear theory alone described by B( 𝑘, 𝜇 ) . A further role of D( 𝑘, 𝜇 ) is to account for any residual broadband term in 𝑃 𝑤 ( 𝑘 ) -ideally this would represent only the BAO wiggles, but in practice this separation cannot be performed on the linear power spectrum exactly (leading to, i.e., differences in the BAO templates using different methods, see §7.4 and Fig. 17). In this case, a good choice of D( 𝑘, 𝜇 ) ensures consistent results regardless of the method used to generate 𝑃 𝑤 ( 𝑘 ) . A particularly simple, well-motivated choice for B( 𝑘, 𝜇 ) is (Kaiser 1987; Peacock & Dodds 1994; Cole et al. 1995)  \nB( 𝑘, 𝜇 ) = ( 𝑏 + 𝑓 𝜇 2 ) 2 GLYPH<18> 1 + 1 2 𝑘 2 𝜇 2 Σ 2 𝑠 GLYPH<19> -2 , (28)  \nwhere in addition to the linear theory prediction we have included a rough 'Fingers of God' (FoG) parametrization through a single free velocity dispersion parameter Σ 𝑠 in order to capture the rapid damping of the power spectrum along the line-of-sight due to smallscale virial motions.  \nPast BAO analyses (e.g., Anderson et al. 2014; Beutler et al. 2017a; Bautista et al. 2021) have tended to rely on parametrizations of D( 𝑘, 𝜇 ) with angular multipoles expanded in polynomials of the wavenumber 𝑘 in powers between some chosen 𝑛 min , max . The logic is that polynomials of sufficiently low order cannot mimic the oscillatory BAO signal. In configuration space, past analyses have also  \nemployed polynomials in 𝑟 of sufficiently low degree with the expectation that the BAO peak in the correlation function cannot be reproduced by polynomials (Ross et al. 2017; Bautista et al. 2021). While historically demonstrated to be robust, this approach suffers from a few limitations: (1) there is no clear criteria for what orders of polynomial terms to include - which is a function of both the scale cuts and noise of the analysis (the polynomial cannot have more critical points than the number of BAO wiggles detected) - and even if there were, there is no clear relation between the order of polynomial and smoothness of the function; and (2) while low-order polynomials cannot mimic wiggles in Fourier space, it is easy for them to reproduce a peak in configuration space. Due to these limitations, the stability of any particular prescription for the broadband given a particular set of data was only empirical. In fact, the earliest broadband models used to fit BAO in the Sloan Digital Sky Survey were based on roughly corresponding powers of 𝑟 𝑛 and 𝑘 -3 -𝑛 to fit the correlation function and power spectrum, respectively (Xu et al. 2012; Anderson et al. 2012) but as the data improved the number of polynomial terms in the power spectrum increased while those in the correlation function stayed roughly constant. In this work we instead develop and validate a new unified treatment of the broadband D( 𝑘, 𝜇 ) for the power spectrum (and correlation function) that avoids these limitations. While the arguments in this section are quite general, we note that their assumptions are implicity and explictly tested in many accompanying DESI 2024 BAO papers-we refer readers interested in the effects of fiducial cosmology and observational systematics on the broadband model developed to Pérez-Fernández et al. (2024); DESI Collaboration et al. (2024a).", '5.1 Fourier Space': "Our starting point will be to develop a broadband model in Fourier space D ℓ ( 𝑘 ) , where we allow for distinct terms for each clustering multipole ℓ . The requirements on D ℓ are that they should not be degenerate with the BAO signal itself; this is really a requirement on the data that any theoretical or observational systematics not have features that can 'look' like the BAO. In practice this requires that D ℓ be smoother than the BAO scale, or rather that its values be uncorrelated beyond some scale Δ > 𝜋 / 𝑟 𝑑 , which is half the wavelength of the BAO by the Nyquist sampling theorem.  \nAconvenient way to do this is to parameterize the broadband using a cubic spline basis with bases separated by a suitably chosen Δ , i.e.  \nD ℓ ( 𝑘 ) = ∞ ∑︁ 𝑛 = -1 𝑎 ℓ,𝑛 𝑊 3 GLYPH<18> 𝑘 Δ -𝑛 GLYPH<19> . (29)  \nHere 𝑊 3 is the 'piecewise cubic spline' (4th-order) extension of the counts-in-cell interpolation kernel (Hockney & Eastwood 1988; Jeong 2010)  \n𝑊 3 ( 𝑥 ) = 1 6         4 -6 | 𝑥 | 2 + 3 | 𝑥 | 3 | 𝑥 | ≤ 1 , 8 -12 | 𝑥 | + 6 | 𝑥 | 2 - | 𝑥 | 3 1 < | 𝑥 | ≤ 2 , 0 | 𝑥 | > 2 , (30)  \nwhich has the convenient property that D ℓ ( 𝑘 > 0 ) = 1 when 𝑎 ℓ,𝑛 = 1. The first few kernels in this basis fit to the linear power spectrum, where we have chosen to use Δ = 0 . 06 ℎ Mpc -1 > 𝜋 / 𝑟 𝑑 , are shown in Figure 6. The model is able to match the broadband of the power spectrum without reproducing the BAO signal itself.  \nA nice feature of this parametrization of the broadband is that the free coefficients 𝑎 ℓ,𝑛 in this model are by construction never degenerate with the BAO signal, and thus cannot lead to numerical instability when fitting the BAO scale. As an extreme example, even  \n0.40  \nFigure 6. Smoothspline basis for broadband fitting. (Top) The basis in Fourier space, multiplied by 𝑎 𝑛 coefficients fit to the linear power spectrum (black solid). The resulting broadband fit, which broadly matches the power spectrum but doesn't contain BAO wiggles, is shown in the black-dashed line. (Middle) The same spline functions without the coefficients 𝑎 𝑛 Hankel-transformed to give the correlation function monopole. The first two basis functions have smooth large-scale contributions while the rest are oscillatory and confined to 𝑟 < 2 𝜋 / Δ . The small wiggles at the largest 𝑟 are due to ringing from the discontinuous higher derivatives of the spline basis. (Bottom) Same as the middle plot but for the quadrupole. The first three basis functions in this case have smooth large-scale contributions.  \n<!-- image -->  \nif 𝑛 Δ > 𝑘 max such that 𝑎 ℓ,𝑛 is unconstrained, it will not interfere with the BAO signal within the range of fit where BAO are detected, but rather simply sample its prior. As such there is no fine-tuning required for the number of free parameters unlike in the polynomial case.  \nFinally, we note that in the Gaussian power spectrum likelihood typically used to constrain the BAO, any linear parameter in the power spectrum model, e.g. the 𝑎 ℓ,𝑛 's, can be analytically marginalized instead of being sampled directly. This is equivalent to adding an additional theoretical contribution to the covariance matrix if one does not care about the actual value of the marginalized-over parameters (Taylor & Kitching 2010, see also Section 6.2). Since we allow  \n𝑎 ℓ,𝑛 to take on arbitrary values, our prescription is equivalent to allowing infinite error at each wavenumber 𝑃 ℓ ( 𝑘 ) , with the proviso that the error across 𝑘 points are correlated when they are within Δ of each other, and the error on different multipoles is uncorrelated. A related theoretical-error based approach to fitting the BAO was proposed in Philcox et al. (2020), who impose an exponentially decaying correlation with Δ of the error but more narrowly estimate its size according to 1-loop perturbation theory. As we will now see, while the spline broadband induces a slightly counterintuitive theoretical covariance that is technically infinite at each 𝑘 point, by performing a simple coordinate (Fourier) transform the theoretical error can instead be strongly localized below some scale in configuration space.", '5.2 Configuration Space': "The spline basis introduced in Fourier space is composed of smooth and well-localized functions which can be readily Hankel transformed into correlation functions. Figure 6 shows their 𝑗 0 , 2 transforms: these basis functions turn into a basis of oscillatory functions in configuration space with zeros around 𝑟 = 𝜋 / Δ , after which they rapidly decay. These functions have arbitrarily high frequencies below this radius as a function of 𝑛 . This is to be expected since the functions in Fourier space have widths given by Δ and are centered at 𝑛 Δ . These functions are given by  \nΔ 3 𝐵 ℓ,𝑛 ( 𝑟 Δ ) = 𝑖 2 ∫ 𝑑𝑘 𝑘 2 2 𝜋 2 𝑊 3 ( 𝑘 Δ -𝑛 ) 𝑗 ℓ ( 𝑘𝑟 ) (31)  \nand we give their exact form in Appendix E. Fitting these functions with free amplitudes in configuration space, i.e.  \n˜ 𝐷 spline ℓ ( 𝑟 ) = Δ 3 ∞ ∑︁ 𝑛 = -1 𝑎 ℓ,𝑛 𝐵 ℓ,𝑛 ( 𝑟 Δ ) , (32)  \nnulls any signal in the correlation function, BAO or otherwise, up to 𝑟 min ≈ 𝜋 / Δ , since the functions sample all available Fourier frequencies up to that scale. While there are, in principle, an infinite number of broadband parameters in this expansion, their effects are localized to small 𝑟 ; we therefore have that the spline broadband model in Fourier space is equivalent to setting an 𝑟 min when performing the fit in configuration space rather than including any smooth broadband (e.g. polynomial) model. The two exceptions to this rule are the two basis functions 𝑛 = [ 0 , 1 ] in the quadrupole, which have some smooth dependence out to large 𝑟 . Note that the 𝐵 2 , -1 term also has nontrivial support at large 𝑘 but it is degenerate with the 𝑘 min contributions discussed below due to being confined to low 𝑘 only. We test the equivalence of these Fourier and configuration space configurations in §7 (see in particular the discussion around Fig. 9). In particular, we note that while there is some 'spillover' due to the finite localization of these functions to small 𝑟 , the equivalence at the level of posteriors is quite insensitive to these numerical choices.  \nPerhaps surprisingly, however, another (implicit) component of the Fourier-space broadband model does make it necessary to introduce a smooth polynomial broadband model in configuration space. In order to protect against unknown large-scale data systematics we typically only include scales above a certain 𝑘 min in the power spectrum when fitting the BAO. The Fourier transforms of these uncontrolled largescale multipoles can be written as an expansion in even powers of the radius 𝑟 (odd powers are cancelled when computing even multipoles), such that  \n˜ D ℓ ( 𝑟 > 𝑟 min ) = ˜ 𝑎 0 + ∞ ∑︁ 𝑛 = 1 ˜ 𝑎 𝑛 GLYPH<18> 𝑟𝑘 min 2 𝜋 GLYPH<19> 2 𝑛 . (33)  \nFor typical 𝑘 min , e.g. 𝑘 min = 0 . 02 ℎ Mpc -1 , then 2 𝜋 / 𝑘 min ≈ 300 ℎ -1 Mpcandthehigherordertermsinthisseries get progressively smaller at the BAO scale. We can hence truncate this to first order, 𝑛 = 1.  \nOverall, both the considerations above lead to a model for the broadband of the monopole and quadrupole of the correlation function given by  \n˜ D 0 ( 𝑟 ) = ˜ 𝑎 0 , 0 + ˜ 𝑎 0 , 1 GLYPH<18> 𝑟𝑘 min 2 𝜋 GLYPH<19> 2 , (34)  \n˜ D 2 ( 𝑟 ) = ˜ 𝑎 2 , 0 + ˜ 𝑎 2 , 1 GLYPH<18> 𝑟𝑘 min 2 𝜋 GLYPH<19> 2 + Δ 3 GLYPH<0> 𝑎 2 , 0 𝐵 2 , 0 ( 𝑟 Δ ) + 𝑎 2 , 1 𝐵 2 , 1 ( 𝑟 Δ ) GLYPH<1>  \n(35)  \nwith 𝐵 2 , 0 ( 𝑟 Δ ) and 𝐵 2 , 1 ( 𝑟 Δ ) given by Eq. 31.", '6 NUMERICAL VALIDATION WITH CONTROL VARIATE SIMULATIONS': 'Having provided a comprehensive theoretical motivation for how the BAO can be modelled and fit given a set of power spectrum or correlation function measurements, we now turn to validating this model using a set of precise simulated measurements with a known cosmological model. In addition to what has been presented so far, there are a number of choices of how exactly one builds the full model that we will also test, many of which are required due to subtle differences between the model we have presented here and what has been done in previous works. The full extent of the things we aim to test/understand are:  \n- · How do the BAO constraints compare between the previous polynomial-based broadband model and the new spline-based model?\n- · What is the impact of different prior choices for the BAO/FoG damping parameters Σ ⊥ , Σ | | and Σ 𝑠 ?\n- · What is the impact of including the FoG term in the wiggle component of the model as well as the smooth component?\n- · What are the optimum scale cuts to use for fitting the power spectrum and correlation function?\n- · How does the choice of method used to perform the wiggle/nowiggle split impact the BAO constraints?\n- · Does dilating the full power spectrum model or just the wiggle component make a difference?\n- · Does our updated model of the correlation function - consistently treating it as the Hankel transform of the dilated power spectrum model with the same bias and RSD parameters; and accounting for the measurement binning - change the BAO constraints compared to previous methods?\n- · Does the choice of MCMC sampler used to perform the fit impact the constraints?  \nCombining all of these effects together, we will arrive at a single consensus value for the modelling systematic error within the DESI BAO analyses, to be combined with the theoretical systematic error from earlier.', '6.1 Simulations': "For all the numerical tests in this work, aimed towards motivating a modelling systematic error budget for DESI, we use a suite of 25 high precision measurements of the power spectrum and correlation function. These are obtained from control variate realisations of a set of cubic simulations populated with a set of galaxies mimicking the  \nclustering of the DESI LRG sample. Here, we give a brief summary of the relevant details of these simulations, and refer the reader to our companion papers, Garcia-Quintero et al. (2024) and MenaFernández et al. (2024) for a more complete description.  \nOur base cubic simulations, the AbacusSummit 8 suite (Maksimova et al. 2021) are a set of over 140 simulations across 97 different cosmological models, each with a box size ( 2 ℎ -1 Gpc ) 3 containing 6192 3 particles; for an effective particle mass of 2 × 10 9 ℎ -1 𝑀 ⊙ . In this work, we restrict to the 'cosm000', Planck Collaboration et al. 2020 Λ CDM-based cosmological model, 9 and a single simulation snapshot with redshift 𝑧 = 0 . 8, which corresponds closely to the peak of the DESI LRG sample number density. The simulations are populated with the 'vanilla' Halo Occupation Distribution model described in Mena-Fernández et al. (2024), leading to a sample with constant number density ∼ 1 × 10 -3 ℎ 3 Mpc -3 . Finally, we apply a reconstruction algorithm to these data, following RecSym convention, and other algorithmic choices presented in (Chen 2024; Paillas et al. 2024).  \nFrom these, a set of more precise clustering measurements is obtained using the control variate (CV) technique (Chartier et al. 2021; Ding et al. 2022; Kokron et al. 2022; DeRose et al. 2023). The CV technique works by subtracting off the intrinsic sample variance in simulated cosmological observables shared with highly correlated but more cheaply computable surrogate observables-in the case of reconstruction, Hadzhiyska et al. (2023) showed that a reasonable such surrogate is the linear 2-point function computed from the same initial conditions. 10 We refer the interested reader to that work for further details, particularly as pertains to the implementation of CV in the context of DESI and AbacusSummit.  \nThe uncertainties on our clustering measurements are obtained using a set of 1000 approximate EZmock simulations (Chuang et al. 2015) with the same volume, particle number density, and HOD as the AbacusSummit simulations. Reconstruction is applied to these approximate simulations in the same way as for the AbacusSummit realisations. Validation of these covariance matrices and comparisons with other approaches for estimating them can be found in Alves & DESI Collaboration (2024); Forero Sánchez & DESI Collaboration (2024); Rashkovetskyi et al. (2024). Because we use a covariance matrix estimated from a finite number of simulations, we also de-bias the inverse covariance matrix and account for noise in its estimation following Hartlap et al. (2007); Percival et al. (2014, 2022).  \nFigure 7 shows the mean power spectrum and correlation function of the 25 realisations along with errorbars for a single realisation and our best fit fiducial BAO model (see §7). This figure shows the extremely precise set of measurements available for us to test our modelling assumptions - and when fitting the mean of 25 realisations, we further reduce the covariance matrix by a factor 25 compared to that shown in the figure to represent the error on the mean. When fitting a single of the 25 individual LRG realisations, the resulting precision on the BAO constraints is ∼ 0 . 3% on 𝛼 iso and ∼ 1% on 𝛼 ap . Fitting the mock mean, with reduced covariance, gives a precision of ∼ 0 . 1% on 𝛼 iso and ∼ 0 . 5% on 𝛼 ap . For comparison, the expected aggregate precision on the BAO constraints for 0 . 0 < 𝑧 < 1 . 1 for DESI Y1 and the full DESI sample are 0 . 51%", '8 https://abacussummit.readthedocs.io/en/latest/': '9 Specifically, with 𝜔 b = 0 . 02237, 𝜔 c = 0 . 12, ℎ = 0 . 6736, 𝐴 s = 2 . 083 × 10 -9 , 𝑛 s = 0 . 9649 and a single massive neutrino with mass 𝑚 𝜈 = 0 . 06 eV. 10 That this is a good choice can be further understood from §4 that shows that the reconstructed density, particularly in RecSym , has its displacement terms reduced to ( 1 - S) ∗ Ψ ( 1 ) . This displacement term is what leads to most of the decorrelation with the initial conditions.  \n(1 . 78%) and 0 . 24% (0 . 68%) for 𝛼 iso ( 𝛼 ap ) respectively (DESI Collaboration et al. 2024d). Hence, a single one of our simulations is sufficient to investigate systematics at the level of ∼ 2 × the precision of DESIs Y1 measurements, while averaging over all 25 realisations is sufficient to determine systematics for the full DESI sample.', '6.2 Numerical implementation in Barry': "To test the various theoretical systematics and modelling choices identified above, we use the modular BAO-fitting code /b.pc/a.pc/r.pc/r.pc/y.pc (Hinton et al. 2020). This code is designed to be fully modular, easily enabling comparative fits to the same data when different modelling choices are made. This is done by specifying a series of unique model and dataset pairs in one script, which are then submitted automatically as unique jobs to a supercomputing cluster. Hinton et al. (2020) used this architecture to study the impact of different perturbationtheory inspired models on the isotropic BAO signal (Seo et al. 2016; Beutler et al. 2017a; Ding et al. 2018; Noda et al. 2020; Chen et al. 2021), finding excellent agreement between the different models. However, as demonstrated in Section 4, we find that a more nuanced and complete treatment of theory leads to a single clear choice of reconstruction convention ( RecSym ), and a single BAO damping term. We hence fix our model to the one presented in Section 5. All of the remaining modelling choices outlined at the start of this Section can be investigated by simple changes to the input parameters within /b.pc/a.pc/r.pc/r.pc/y.pc.  \nTo facilitate this, significant improvements have been made to the code relative to that presented in Hinton et al. (2020) which we will now describe. For all other base aspects of the code, we refer the reader to this reference. In contrast to Hinton et al. (2020), the latest version of /b.pc/a.pc/r.pc/r.pc/y.pc implements anisotropic BAO fitting as default, reverting back to a standard isotropic method only when specified for particular model-dataset pairs.  \nIt is important to note that an 'isotropic' BAO fit, both within /b.pc/a.pc/r.pc/r.pc/y.pc and in general, is not the same as running an 'anisotropic' fit to only the monopole with fixed 𝜖 = 0. The reason for this is due to the impact of the survey window function. Both the window function convolution and modelling of wide-angle effects are implemented via matrix multiplications as described in Beutler & McDonald (2021), which also allows the fitting of the odd-multipoles of the power spectrum if they are provided. Even when fitting only the monopole, the presence of a survey window function means that power leaks from higher-order multipoles into the monopole and vice versa, requiring the full set of multipoles to be modelled and the window function for all even and odd multipoles to be used. When an isotropic fit is performed, the higher-order multipoles are assumed to be zero and the part of the window function responsible for describing the leakage of power between these and the monopole is ignored. This is a fundamental difference in the modelling that is retained in /b.pc/a.pc/r.pc/r.pc/y.pc for historical reasons and comparison to results from previous surveys, but in this work, we only consider the correct 'anisotropic' approach of modelling all higher-order multipoles regardless of the exact parameters, scales or data multipoles being fit.  \nA second feature that is important to highlight is the ability to analytically marginalise over the broadband terms during the BAO fitting. As these enter at linear order, such marginalisation can be done fully (by modifying the likelihood used during MCMC fitting or optimisation) or partially (by using a standard Gaussian likelihood and finding the best-fit broadband terms at each step of the fitting). Appendix F discusses this further and demonstrates the consistency between these approaches. All results presented hereon use full analytic marginalisation, which we argue is the formally correct choice,  \nFigure 7. The mean clustering signal, taken from 25 control variate simulations of DESI's LRG sample, used for all the modelling tests in this work. In all panels closed/open points show the mean measured monopole and quadrupole clustering respectively, and error bars are those for a single realisation (i.e., not the error on the mean). Left panels are the power spectrum and right are the correlation function. For both, the top panels show the full signal, including the best-fitting full- and no-wiggle models (solid/dashed lines respectively) using our fiducial fitting methodology (see Section 7). The bottom panels show only the remaining wiggle component of the multipole data/models after subtracting the no-wiggle model (i.e., showing only the second term on the right-hand side of Eq. 5), multiplied by 𝑘 , and arbitrarily offset by ± 100 ( ± 10) for the power spectrum (correlation function) for visualisation purposes.  \n<!-- image -->  \nalthough partial analytic marginalisation is retained as an option in the code for comparison to previous works which used this approach (i.e., Ross et al. 2017; Bautista et al. 2021).  \nThe full code used to perform all the tests in this work is available at https://github.com/Samreay/Barry .", '7 RESULTS': 'This section presents the results of our various systematic modelling tests. Unless otherwise stated, we fit the post-reconstruction data with Gaussian priors on the BAO and FoG damping parameters, Σ ⊥ = N( 2 , 1 ) , Σ | | = N( 5 , 2 ) , Σ 𝑠 = N( 2 , 2 ) on scales 0 . 02 ℎ Mpc -1 ≤ 𝑘 ≤ 0 . 30 ℎ Mpc -1 and 50 ℎ -1 Mpc < 𝑠 < 150 ℎ -1 Mpc for the power spectrum and correlation function respectively. Both the linear galaxy bias and RSD parameters are allowed to vary within wide flat prior ranges. We dilate only the wiggle component of the model, and the FoG damping applies only to the smooth component. The wiggle/no-wiggle split is performed using the method of Wallisch (2018) and all our fits use the Nautilus sampler (Lange 2023).', '7.1 Validating the spline-based broadband model': 'Our first set of tests validates the new spline-based broadband method. First, we test the dependence of our fit results to the smoothing scale Δ . Intuitively, splines with too small a Δ will be able to reproduce the BAO signal, while splines with too large a Δ will lack sufficient freedom to fit the broadband. This is especially clear in configuration space where the former corresponds with 𝑟 min overlapping the BAO peak, while the latter corresponds with 𝑟 min extending too  \nfar out from the BAO peak. Fortunately, no fine tuning in Δ is required as fits using the spline method are very stable over a range of Δ until the BAO scale is reached.  \nFigure 8 shows fits to the mean LRG power spectrum with varying values of Δ . The constraints on 𝛼 ∥ , ⊥ are extremely stable to a broad range of Δ , and only begin to rapidly degrade around Δ ∼ 0 . 035 ℎ Mpc -1 or roughly 10% above the Nyquist frequency of the BAO wiggles. The central values of the fits shift by less than 0 . 1% in all cases and the best-fit 𝜒 2 changes only by around the amount implied by the increase in free spline parameters with support within the fitting range ( 𝑘 max = 0 . 3 ℎ Mpc -1 ). Since there is no apparent bias due to using a slightly wider spacing Δ but a clear degradation of the overall signal well before the Nyquist frequency is reached we pick roughly twice this frequency, Δ = 0 . 06 ℎ Mpc -1 , as our fiducial choice. At 𝑘 max = 0 . 3 ℎ Mpc -1 this corresponds to 7 free parameters per multipole.  \nNext, we confirm that the Fourier-space broadband model is equivalent to setting 𝑟 min = 𝜋 / Δ with a minimal set of spline coefficients and large-scale polynomial terms. To do so we perform two fits: one where the spline broadband including 20 bases with Δ = 0 . 06 ℎ Mpc -1 is Hankel transformed and the correlation function fit from 0 ℎ -1 Mpc < 𝑟 < 130 ℎ -1 Mpc (summing from 𝑛 = -1 to 𝑛 = 19 in Eq. 32), and another where only the first two spline bases in the quadrupole with large-scale support in 𝜉 2 are kept (using only the 𝑛 = 0 and 𝑛 = 1 terms in Eq. 32), fitting to 80 ℎ -1 Mpc < 𝑟 < 130 ℎ -1 Mpc, where 80 ℎ -1 Mpc is chosen as the 𝑟 min where the spline basis has essentially no support. The two constraints are shown in Figure 9. The BAO constraints are essentially identical, up to percent level differences in constraining power, despite the radial fitting range differing by close to a factor of two. Both the means and variances of the posteriors on the linear theory amplitudes 𝑏 and 𝑓 also agree at the roughly 5% level, which is  \nFigure 8. Fits to the mean power spectrum with varying levels of Δ . The constraints on the BAO parameters 𝛼 ∥ , ⊥ are remarkably stable to Δ up to 𝑘 ∼ 0 . 04 ℎ Mpc -1 ; at around Δ ∼ 0 . 035 ℎ Mpc -1 , just shy of 10% above the Nyquist frequency, the constraints rapidly degrade.  \n<!-- image -->  \nnotable given the very different radial ranges over which the data are fit. This suggests that these linear-theory amplitude parameters are mostly fitting the amplitude of the BAO only while sensitivity to the amplitude of the smooth component is removed either by 𝑟 min or by fitting the spline basis.  \nFinally, we would like to check that the fits are robust against various levels of noise. In order to perform this check we fit the data with the covariance re-scaled by factors of 𝐴 = 1 / 4 , 1 / 2 , 1 , 2 , 4. This corresponds approximately to rescaling the volume by the inverse factor - at 𝐴 = 1 / 4 and 𝐴 = 4 these rescalings correspond to more volume than all the DESI Y1 samples combined and less volume than BOSS CMASS. The results in Fourier and configuration space are shown in Figure 10. The fits are remarkably stable across these rescalings, with the error bars on 𝛼 ∥ , ⊥ scaling as √ 𝐴 to a very good approximation, demonstrating that the spline broadband does not sap BAOinformation out of the signal at any SNR level. Indeed, we have tested that this remains true even for a 16 × rescaling, beyond which the BAO scale can essentially only be constrained to order unity.', '7.2 Impact of BAO/FoG damping parameter choices': "The BAO damping parameters Σ | | and Σ ⊥ play a key role in describing the shape of the measured BAO feature, which in turn impacts how well its peak position/wavelength can be determined. Furthermore, Σ 𝑠 is designed to account for non-linear RSD in the broadband shape of the clustering but may interact with the BAO signal depending on how it is implemented. Previous studies (Alam et al. 2017; Gil-Marín et al. 2020; Bautista et al. 2021) have typically fixed the values of these parameters based on fits to a set of mocks, additionally demonstrating that fits to the data remain robust if different fixed values are used, or the parameters are freed up entirely. However, the treatment of Σ 𝑠 has not always been consistent, and it is not clear to what extent the priors on the three Σ parameters impact BAO constraints at DESI precision.  \nTo test, this we first investigate how fixing these to different  \nFigure 9. Fits to the mean correlation function using a wide scale range 0 ℎ -1 Mpc < 𝑟 < 130 ℎ -1 Mpc and the full set of Hankel-transformed spline coefficients compared to when an 𝑟 min cut of 80 ℎ -1 Mpc is applied and only the 𝑛 = [ 0 , 1 ] spline bases in the quadruopole are used. The excellent consistency between the two demonstrates that our default method of fitting the BAO removing small scales and unnecessary spline components is sufficient.  \n<!-- image -->  \nvalues changes the BAO constraints. Figure 11 shows the bias Δ 𝛼 iso = 𝛼 iso -1 and similar for 𝛼 ap when we fit the mock mean clustering. Our results show that fixing the choice of damping parameters is somewhat risky. While the constraints on the BAO dilation parameters are remarkably robust to moderate changes in the choices of damping scales, both the power spectrum and correlation function show systematic variation of increasing bias in the 𝛼 's when the BAOdamping Σ 's are increased significantly - biases of order 0 . 1% and 0 . 2% for 𝛼 iso and 𝛼 ap respectively occur when the Σ values are changed by more than ∼ 2 ℎ -1 Mpc from their fiducial values. Similar trends have been identified in previous works, see e.g., Table B2 in Ross et al. (2017). We see no notable impact from varying Σ 𝑠 on either the correlation function or power spectrum, as expected given our default model decouples it from the wiggle power spectrum.  \nThe bias is more pronounced for the combination of 𝛼 iso with Σ ⊥ than with Σ | | , while 𝛼 ap shows a complex behaviour from variations in either damping parameter. This can be understood from theoretical considerations - Taylor expanding the full BAO damping term in Eq. 26 and integrating with respect to the ℓ = 0 and ℓ = 2 Legendre polynomials allows us to write effective BAO damping scales for the monopole and quadrupole as  \nΣ 2 ℓ = 0 = 1 1 + 2 𝛽 + 3 5 𝛽 2 GLYPH<20> GLYPH<18> 1 + 6 5 𝛽 + 15 35 𝛽 2 GLYPH<19> Σ 2 | | + GLYPH<18> 2 + 4 5 𝛽 + 6 35 𝛽 2 GLYPH<19> Σ 2 ⊥ GLYPH<21> , (36) Σ 2 ℓ = 2 = 1 2 𝛽 + 6 7 𝛽 2 GLYPH<20> GLYPH<18> 1 + 12 7 𝛽 -5 7 𝛽 2 GLYPH<19> Σ 2 | | -GLYPH<18> 1 -2 7 𝛽 -1 7 𝛽 2 GLYPH<19> Σ 2 ⊥ GLYPH<21> , (37)  \nwhere 𝛽 = 𝑓 / 𝑏 is the linear RSD parameter. Information on 𝛼 iso arises primarily from the monopole and so bias in this parameter actually occurs when variations in Σ | | or Σ ⊥ cause a change in Σ ℓ = 0 . For the LRG mocks used here, where 𝛽 ≈ 0 . 35, the same constant  \nFigure 10. Fits to the mean power spectrum (top) and correlation function (bottom) when the covariance matrix is re-scaled between 0.25x and 4x. The excellent consistency of the constraints demonstrates that our default splinebasis never removes BAO information regardless of the SNR level.  \n<!-- image -->  \nchange applied to Σ 2 ⊥ has a ∼ 1 . 5 × greater impact on the damping scale in the monopole, compared to the same change applied to Σ 2 | | . This is consistent with the larger bias seen in 𝛼 iso when Σ ⊥ is varied compared to Σ | | . Conversely, the damping scale in the quadrupole is more affected by changes in Σ | | , but can also become too small if Σ ⊥ is increased enough at fixed Σ | | . As the quadrupole is the primary source of information on 𝛼 ap , considering variations in Σ 2 ℓ = 2 can explain the more nuanced behaviour of the 𝛼 ap as a function of both Σ | | , ⊥ . We note that the dependence on 𝛽 meansthat the picture will be slightly quantitatively different for other tracers, though qualitatively similar.  \nOverall, this suggests that some freedom should be given to the BAO damping parameters. We recommend using a Gaussian prior, centred on theory calculations such as in Fig. 5 or from fitting sim-  \nFigure 11. Percentage bias in the BAO dilation parameters as a function of the fixed damping scale used in the model when fitting the mean of the 25 LRG simulations. Each band shows the mean and standard deviation in the BAO parameters (relative to the expected value of 1) from fitting our mock mean measurements. Left panels show fits to the correlation function, while right panels are fits to the power spectrum. In all cases, the vertical lines show our fiducial choice of damping scales, while the horizontal lines are placed at characteristic values of 0 . 1% for 𝛼 iso and 0 . 2% for 𝛼 ap to guide the eye. These results suggest that the BAO model can accommodate small inaccuracies in the fiducial choice when fixing the damping scales, but that larger-than-ideal biases occur if the fiducial value is more than ∼ 2 ℎ -1 Mpc from the correct choice. This motivates allowing for some freedom in the damping scale, i.e., a Gaussian prior centred on the fiducial value with width 1 -2 ℎ -1 Mpc.  \n<!-- image -->  \nulations, but with width 1 -2 ℎ -1 Mpc to allow for the presence of non-linearities not captured by theory or deviations between the mocks and data which may otherwise cause small 0 . 1 -0 . 2%biases. When fitting the damping scale to simulations we advocate for fitting the BAO in high-precision simulated mocks (e.g. control-variate DESI mocks) using Equation 26 with the cosmology set to 'truth'  \n(i.e. 𝛼 ∥ , ⊥ = 1) - previous works often instead fit to the 'propagator', i.e. the correlation between the initial conditions and (reconstructed) galaxy density, but we note that this quantity is not strictly the same as the BAO damping since it is not cutoff at the BAO scale as in Equation 9.  \nAn alternative is to allow the BAO damping parameters to freely vary between wide flat priors such that the values are purely informed by the data. However, in other tests as part of this work, on more weakly constraining simulations than the control variates used here, we found that this generally reduces constraining power and has the potential to bias the BAO constraints. In the literature, nearly all previous BAO studies have elected to fix the damping parameters; of the studies we referred to only Gil-Marín et al. (2020) tested varying these, finding consistent BAO constraints, but we caution that the combined BOSS+eBOSS LRG dataset used there is particularly well constraining. We also note that future studies could also investigate whether it is more useful to place priors on Σ ℓ = 0 and Σ ℓ = 2 rather than Σ ⊥ and Σ | | as our results suggest that it is changing along these degeneracy directions that ultimately cause biases in the BAO scale.  \nHaving determined a reasonable choice of prior for the damping parameters, we then test for any impact of including the FoG damping term in the wiggle component of the model as well as the smooth component, i.e., including the Lorentzian term from Eq. 28 in Eq. 26. As a first demonstration of why we advocate for decoupling the FoG damping from the wiggle component of the model, we show in Fig. 12 that doing so effectively removes the degeneracy between Σ | | and Σ 𝑠 when fitting the mean of our mocks. This makes it far easier to enforce or determine suitable priors on these parameters in the manner described above and improves the convergence of the fitting. This breaking of the degeneracy occurs even when the 𝛼 's are allowed to vary in concert with the damping parameters.  \nSecondly, we look at the impact of this on the BAO constraints themselves by fitting each of our 25 individual mock realisations and then looking at the distribution of mock-to-mock differences in the 𝛼 's for the two different modelling configurations. The results are shown in Fig. 13, where we find that regardless of the choice of clustering statistic, or whether a polynomial-based or our new spline-based broadband model is used, the largest average difference is ≤ 0 . 02% and 0 . 11% for 𝛼 iso and 𝛼 ap respectively. Albeit small, this latter difference is detected with 2 𝜎 significance and so should be accounted for in our final systematic error budget. That said, the largest difference seen in any individual realisation is ∼ 0 . 25% (0 . 75 𝜎 ) for 𝛼 iso and 0 . 75% (0 . 75 𝜎 ) for 𝛼 ap , indicating that while on average the results are robust to the choice of where the FoG damping is included, this choice could lead to statistical fluctuations on individual datasets dependent on the noise in the data.", '7.3 Impact of different scale cuts': 'We next turn to verify the appropriate range of scales to include in the BAO fit. For both the power spectrum and correlation function we fit the mock mean clustering with a variety of choices for the minimum and maximum scale cuts and measure the bias in the dilation parameters, and their uncertainties relative to our default. These are shown in Figures 14 and 15 for the power spectrum and correlation function respectively.  \nFor the power spectrum, we find that the difference from unity exceeds 0 . 1% and 0 . 2% for 𝛼 iso and 𝛼 ap respectively for a couple of narrow ranges of 𝑘 max . The 𝑘 max at which this occurs is not the same for the two dilation parameters, and, more importantly, including even more non-linear scales does not make the bias worse. There is no significant trend of increasing bias when increasing/decreasing  \nFigure 12. Constraints on the BAO and FoG damping parameters when fitting the mock mean correlation function with fixed 𝛼 iso = 𝛼 ap = 1 and including the FoG damping in only the smooth model component ( 𝜉 ( 𝑠 ) FoG Smooth), or allowing it to multiply both the smooth and wiggle components ( 𝜉 ( 𝑠 ) FoG Wiggles). The dashed lines represent the central values we use for the priors on these parameters in our other fits. Including the FoG damping in only the smooth model component is preferable as it removes a clear degeneracy between the BAO and FoG damping scales.  \n<!-- image -->  \neither 𝑘 min or 𝑘 max . As such, we attribute these regions with slight deviations from unity to noise in the simulations, rather than a bias in the model - looking at the measurements in Fig. 7, one can see a slight excess of power at around 𝑘 = 0 . 29 ℎ Mpc -1 which inconveniently coincides with our fiducial choice of scales.  \nTurning to the uncertainty on the BAO parameters using the power spectrum, we do see a clear weakening of the constraints by up to 50% when restricting the fitting to larger scales 𝑘 max = 0 . 2 ℎ Mpc -1 and some small 10% improvement when including smaller scales 𝑘 max = 0 . 4 ℎ Mpc -1 . There is also perhaps a very weak trend as 𝑘 min is increased and starts to encroach on the first BAO wiggle. It is worth reminding that as this is from a fit to the mean mock measurements, the uncertainties on the clustering measurements are significantly smaller than we expect for DESI DR1 and comparable to the aggregated BAO across 0 . 0 < 𝑧 < 1 . 1 for the full DESI survey. This exacerbates the trends in the uncertainties, and so we conclude that our fiducial choice of fitting scales 0 . 02 ℎ Mpc -1 < 𝑘 < 0 . 30 ℎ Mpc -1 is reasonable - it is not worth extending the fit to smaller scales for DESI DR1 or its upcoming DESI 2024 results. But, as we see no evidence of additional bias, this may be worth testing for the full DESI data.  \nFor the correlation function, we find that the results are remarkably stable to the choice of scale cuts, with a small increase in bias but a decrease in error as 𝑠 min is reduced. The bias reaches 0 . 1% in 𝛼 iso for 𝑠 min < 28 ℎ -1 Mpcacross the range of 𝑠 max we test, and ∼ 0 . 15% in 𝛼 ap when a maximum fitting scale 𝑠 max < 130 ℎ -1 Mpc is also applied. When the upper fitting range is reduced to this, there is enough freedom in the quadrupole model (particularly in its amplitude/BAO damping) that one starts to affect the BAO information if additional scales are not included on either side of the BAO feature with which  \nFigure 13. Difference in the BAO constraints between our 25 mock realisations fit with and without the FoG damping applied to the wiggle component in addition to the smooth component (without is our default). Each box-and-whisker plot shows the median difference over the individual mocks (green line) the interquartile range (blue/orange boxes for power spectrum/correlation function respectively) and the largest individual differences (black error bars). These plots are shown for both the polynomial and spline-based broadband methods. The top two panels show the percentage differences between the same realisation, while the bottom two panels show the difference relative to the average statistical uncertainty on the BAO parameter for a single realisation. The horizontal lines are to guide the eye and are placed at 0.1 for Δ 𝛼 iso and the error-normalised differences, and 0 . 2%for Δ 𝛼 ap .  \n<!-- image -->  \nto constrain the broadband model. Beyond this, any choice of fitting scale seems as good as any other, with < 10% relative change in the uncertainties across all choices, so we opt to fix to our fiducial choice of 50 ℎ -1 Mpc < 𝑠 < 150 ℎ -1 Mpc.', '7.4 Impact of de-wiggling methodology': "When performing a BAO fit, a choice must be made of how to isolate the wiggle and no-wiggle components in the linear power spectrum. A number of different algorithms/approaches for this have been developed; using fitting formulae (Eisenstein & Hu 1998), filters designed for signal processing (such as the Savitsky-Golay filter, B-splines, or a Gaussian filter as used in Hinton et al. 2017), or methods that numerically search for changes in the gradient in the power spectrum or its Fourier transform (Wallisch 2018; Brieden et al. 2022). As shown in Fig. 16, the wiggle component of the  \nFigure 14. The impact of varying minimum and maximum scale cuts ( 𝑘 min and 𝑘 max ) on BAO dilation parameters extracted from the mock mean LRG power spectrum. The top panel is coloured according to the percentage deviation of the fitted 𝛼 's from unity. The bottom panel is the uncertainty of the 𝛼 's relative to our fiducial fitting range (indicated by the green cross). In both cases, the left (right) panel is for 𝛼 iso ( 𝛼 ap ). The regions enclosed by the dashed line in both of the left panels are where the deviation in 𝛼 iso from unity is < 0 . 1%. For the right panels, it denotes a deviation in 𝛼 ap of < 0 . 2%.  \n<!-- image -->  \nmodel arising from these different choices of algorithm can differ appreciably and, without a clear theoretical reason to prefer one over the other, the impact of this choice must be explored.  \nTo do this, we follow the same procedure used previously and fit our 25 individual LRG mock realisations using each of the different template extraction methods. We then take the mock-to-mock difference between the constraints relative to our fiducial Wallisch (2018) method. The results are shown in Fig 17. Despite the clear differences in the templates, the systematic difference in the BAO constraints is < 0 . 02% ( < 0 . 09%) for 𝛼 iso ( 𝛼 ap ). However, as with the case where we investigate the impact of including the FoG damping in both the smooth and wiggle model components, individual realisations can show larger statistical variations when the template is varied - up to 0 . 2% (0 . 75 𝜎 ) for 𝛼 iso , and 0 . 5% (0 . 5 𝜎 ) for 𝛼 ap . The correlation function appears slightly more sensitive to these variations than the power spectrum, consistent with the larger variations between templates over our fitting range seen in the correlation function of Fig. 16.  \nFigure 15. As for Fig. 14, but showing BAO constraints from the correlation function when varying the minimum and maximum fitting scales 𝑠 min and 𝑠 max .  \n<!-- image -->  \nFigure 16. A comparison of the BAO template extracted using a variety of different literature approaches applied to the same linear power spectrum.  \n<!-- image -->", '7.5 Impact of dilating the smooth model component or not': 'BAOfitting constraints arise from dilating the scales in the two-point clustering signal by the 𝛼 parameters and then comparing the dilated model (Eq.5) to the data. In our fiducial set-up, we dilate only the wiggle term 𝐶 ( 𝑘, 𝜇 ) 𝑃 𝑤 ( 𝑘 ) , evaluating the remainder of the model at the undilated scales/angles 𝑘 and 𝜇 . This ensures that constraints', 'Template Dewiggling Differences': "Figure 17. As for Figure 13, but showing the differences in the BAO constraints between our 25 mock realisations using different BAO template extraction methods. In all cases, the template cosmology and linear power spectrum are the same and the differences are taken relative to our fiducial (Wallisch 2018) method.  \n<!-- image -->  \non the 𝛼 's arise only from the wiggle-component. However, one could also apply the dilation to the smooth component of the model. Previous results in the literature have often taken this latter approach (e.g., Beutler et al. 2017a; Ross et al. 2017; Gil-Marín et al. 2020; Bautista et al. 2021).  \nFollowing the same approach as above, fitting the 25 individual LRG mocks with both approaches and looking at the distribution of differences in the BAO constraints, we find that the choice of whether or not to dilate the smooth model component in addition to the wiggles results in a systematic error of < 0 . 02% ( < 0 . 09%) for 𝛼 iso ( 𝛼 ap ). The largest individual differences are ∼ 0 . 15% (0 . 5 𝜎 ) for 𝛼 iso , and 0 . 75% (0 . 75 𝜎 ) for 𝛼 ap .", '7.6 Impact of different methods for transforming from the power spectrum to the correlation function': 'In this work, we provide a fully consistent/equivalent model for the correlation function and power spectrum, treating the former purely as the Hankel transform of the latter, and applying all our modelling choices and BAO dilation to the power spectrum. The new spline broadband method also ensures a consistent modelling between the two. However, some previous analyses (Anderson et al. 2014; Ross et al. 2017) instead transformed the undilated power spectra, without  \ngalaxy bias, and applied both the BAO dilation and galaxy bias terms in configuration space. They also used two different galaxy bias terms for the two correlation function multipoles rather than including a physical angular-dependent RSD with free parameter 𝑓 or 𝛽 . Eqs. 1317 and accompanying text in Ross et al. (2017) describe this process.  \nWeestimate the impact of the differences between these two methods by implementing both and fitting the individual LRG realisations using the spline-based broadband. We find a negligible systematic difference in 𝛼 iso , but a detection of a small difference of 0 . 12 ± 0 . 03% in 𝛼 ap . The fact that the two approaches differ mainly in their modelling of the quadrupole and the difference occurs in 𝛼 ap indicates the presence of a real difference, not just noise. Although we believe that our new approach is preferable, we hence accommodate for this difference in our total systematic error on 𝛼 ap . On individual realisations, we find that the largest difference between the two methods is 0 . 04% (0 . 15 𝜎 ) for 𝛼 iso , and 0 . 4% (0 . 4 𝜎 ) for 𝛼 ap . Despite the small detection of systematic differences when averaging over all 25 realisations, this is consistent with the statistical differences on individual realisations seen from other modelling choices.', '7.7 Impact of different parameter sampling methods': "The last choice we investigate is less about the evaluation of the modelling for the BAO, but rather how the constraints on the 𝛼 's are obtained numerically. /b.pc/a.pc/r.pc/r.pc/y.pc exposes a number of well-known sampling algorithms with which one can fit the data and recover posterior samples (and possibly evidences) on the dilation parameters - these include the publicly available packages /e.pc/m.pc/c.pc/e.pc/e.pc, /d.pc/y.pc/n.pc/e.pc/s.pc/t.pc/y.pc, zeus and /n.pc/a.pc/u.pc/t.pc/i.pc/l.pc/u.pc/s.pc (Foreman-Mackey et al. 2013; Speagle 2020; Karamanis et al. 2021; Lange 2023), each with automated convergence checking.  \nFig. 18 shows the marginalised posterior means on the 𝛼 's for our individual realisations for each of the samplers we test, relative to those obtained using /n.pc/a.pc/u.pc/t.pc/i.pc/l.pc/u.pc/s.pc. We also investigate the width of the 68% and 95% confidence intervals on the 𝛼 's. Given these samplers are all well tested and use reasonable automated convergence checks we do not expect and do not see any systematic differences. What is interesting to quantify however is that there is statistical differences between the posterior means for individual realisations depending on the choice of sampler - up to 0 . 1% and 0 . 5% on 𝛼 iso and 𝛼 ap respectively. In all cases, we see much better convergence in the 68% and 95% confidence intervals than the posterior mean, consistent with statistical sampling noise around the peak of the posterior.  \nWe measured the speed to convergence of these various samplers across the 25 realisations and found that /n.pc/a.pc/u.pc/t.pc/i.pc/l.pc/u.pc/s.pc performed significantly faster than other methods. It also has the benefit of returning the Bayesian evidence as well as posterior samples. For these reasons, we adopted /n.pc/a.pc/u.pc/t.pc/i.pc/l.pc/u.pc/s.pc as our default method, including in this paper.", '7.8 Summary of methodology systematics and differences between our new method and the (e)BOSS methodology': "Overall, our tests in this Section have shown that BAO measurements are remarkably robust to model choices not specified by theory. We find that, if our suggested best-practices (summarised further below) are followed, we expect a systematic error of no more than 0 . 05% in 𝛼 iso and 0 . 12% in 𝛼 ap to arise from any modelling choices. For clarity compared to past studies and as a guide for the future, we now reiterate our best-practices, highlighting differences with previous literature:  \nFigure 18. Difference in the BAO constraints between our 25 mock realisations fit with different MCMC/Nested samplers, relative to those obtained using /n.pc/a.pc/u.pc/t.pc/i.pc/l.pc/u.pc/s.pc. Each box-and-whisker plot shows the median difference over the individual mocks (green line) the interquartile range (blue/orange boxes for power spectrum/correlation function respectively) and the largest individual differences (black error bars). The top two panels show the percentage differences between the same realisation, the middle two panels show the differences between the width of the 68% confidence intervals (i.e., in the 1 𝜎 bounds), and the bottom two panels show the same for the 95% confidence intervals. The horizontal lines are to guide the eye and are placed at 0.1 for Δ 𝛼 iso and 0 . 2% for Δ 𝛼 ap .  \n<!-- image -->  \n(i) For BAO reconstruction, previous works (Xu et al. 2012; Anderson et al. 2012, 2014; Alam et al. 2017, 2021) largely followed RecIso convention. In contrast, we have demonstrated a clear theoretical preference for the RecSym reconstruction convention (see Section 4.1 for a description of these two conventions, and Eqs. 24 and 26 for the corresponding RecSym propgators). Our companion papers Chen (2024) and Paillas et al. (2024) further detail how this can be applied optimally in the context of DESI.  \n(ii) We use the spline-based broadband model validated in §7.1. Previous analyses used polynomials (of slightly varying degrees and functional forms) to fit the broadband. We find the spline method has better theoretical motivation, and to be less arbitrary in its choice of free parameters. For our default fitting of the monopole and quadrupole, between the ranges of 𝑘 = 0 . 02 -0 . 30 ℎ Mpc -1 and 𝑠 = 50 -150 ℎ -1 Mpc and with a spline-width Δ = 0 . 06 ℎ Mpc -1 , this gives 14 free broadband parameters for the power spectrum and 6 for the correlation function. We again clarify that the broadband between the two has the same theoretical basis - the difference in  \nnumber of parameters stems only from the behaviour of the Hankel transform and the different fitting ranges used for the two sets of clustering measurements.  \n(iii) WerecommendchoosingtheBAOdampingparametersbased on theoretical calculations and allowing these to vary within reasonably tight priors (1 -2 ℎ -1 Mpc).ForDESI,wegenerallyuseGaussian priors of width 1 ℎ -1 Mpc for Σ ⊥ and 2 ℎ -1 Mpc for Σ | | and Σ 𝑠 , corresponding to approximately 20% variations from their expected values. Central values for these priors are chosen based on fits to simulated data and theoretical considerations. More info on the priors adopted for our fiducial data analysis can be found in DESI Collaboration et al. (2024a). We found that completely fixing the BAO damping has the potential to lead to biased 𝛼 's if the wrong values are chosen, while allowing them to vary freely can also cause bias, and generally weak constraints, in the case of noisy data. Accounting for shot-noise and non-linearities is able to give robust analytic predictions. Previous BOSS/eBOSS analyses fixed these based on fits to simulations (although found consistent results compared to when these were allowed to freely vary).  \n(iv) We apply the Finger-of-God damping to only the smooth model components, allowing for no interaction with the BAO wiggles. This makes interpreting the impact of the BAO damping parameters simpler. Previous models applied the FoG term equally to both the wiggle and no-wiggle components.  \n(v) The DESI BAO template is constructed following the method of Wallisch (2018) Previous BOSS and eBOSS analyses have used varying methods for extracting the template (i.e., Eisenstein & Hu 1998; Kirkby et al. 2013). Our tests do not point to a single correct choice, but we find a maximum systematic error associated with this choice of only 0 . 02% and 0 . 09% for 𝛼 iso and 𝛼 ap respectively when following our fiducial fitting methodology.  \n(vi) We dilate only the BAO wiggles by the 𝛼 's. Previous analyses generally dilated both the smooth and wiggle components of the model. We believe the new method to be preferable as it guarantees only information from the BAO is captured. Given the large amount of freedom in the broadband used in previous analyses, we do not expect any information was included from the full-shape of the clustering with the previous method, but our model choice ensures this remains true even if the broadband model is less flexible.  \n(vii) We give both the galaxy bias 𝑏 and RSD term 𝛽 flat priors, finding that difficulties determining the appropriate central value of other prior choices can lead to poor fits on (a small number) of individual realisations. Previous analyses typically used more informative (but reasonably wide) Gaussian priors on one or both of these.  \n(viii) We provide a fully consistent/equivalent model for the correlation function and power spectrum, treating the former purely as the Hankel transform of the latter, and applying all our modelling choices and BAO dilation to the power spectrum. BOSS analyses instead transformed the undilated power spectra, without galaxy bias, and applied both the BAO dilation and galaxy bias terms in configuration space. They also used two different galaxy bias terms for the correlation function multipoles rather than the RSD parameter 𝛽 . We find the new, consistent, modelling method to be preferable, and find 0 . 12% of systematic difference in 𝛼 ap between the two methods.", '8 CONCLUSIONS': "Baryon acoustic oscillations (BAO) are one of the premier probes of the cosmic expansion history, especially at high redshifts, providing percent-level and below measurements of cosmological distances up  \nto 𝑧 = 2. Since BAO are an oscillatory signal in Fourier space, with a distinct and fast oscillation frequency set by the sound horizon 𝑟 𝑑 at the baryon drag epoch, its imprint on large-scale structure cannot easily be reproduced by either survey systematics or nonlinearities, making the measured BAO scale remarkably robust. Indeed, in the galaxy power spectrum, the latter typically induces shifts of no more than a percent, shifts which are further reduced by reconstruction. On the other hand, modern spectroscopic surveys, both planned and currently online, promise to probe galaxy clustering with larger volumes and greater number densities than ever before, leading to BAO measurements with statistical power verging upon the bias naively induced by these systematic effects. Indeed, the Dark Energy Spectroscopic Instrument (DESI), which will soon be releasing its first data release (DR1) is projected to yield a cumulative 0 . 2% measurement of the BAO scale when its five-year run is completed.  \nGiven the above, our goal in this paper has been to perform a thorough accounting of the theory and modeling systematics incurred in measuring the nonlinear BAOingalaxy clustering. Our approach has been two-fold: first, we reviewed the theoretical literature modeling the BAO signal in galaxies, filling in holes and updating results where necessary, in order to both estimate the shift in the BAO scale due to nonlinearities and suggest an appropriate fitting form for the observed signal. Secondly, we test these methods on a set of high-precision mocks whose sample-variance has been reduced using Zeldovich control variates, paying particular attention to systematic effects of numerical choices not uniquely specified by the theory.  \nIn the first part of the paper we (re)examine the perturbation theory (PT) of the BAO, particularly as pertains to the nonlinear shift of the nonlinear BAO peak and its damping by long-wavelength displacements. Two of the largest contributions to the former are (a) a nonlinear 'contraction' of the BAO due to density modes larger than 𝑟 𝑑 , whose effect in redshift space we derive and (b) imprints of relative baryon-dark matter perturbations on galaxy clustering. The former effect is substantially larger in redshift space, where spectroscopic surveys naturally operate, than previous calculations in real-space suggest, but is significantly diminished after reconstruction, up to an error induced by a mismatch of the fiducial and true matter clustering amplitude Δ 𝜎 2 8 / 𝜎 2 8 . The impact of these two nonlinear effects can computed as a function of bias and redshift and added to the theoretical error budget of BAO measurements. Post reconstruction, assuming that the matter clustering on large scales ∝ 𝜎 8 is known to 5%, effect (a) leads to shifts in the BAO scale of no more than 0 . 1% in the regimes relevant for DESI, while effect (b) leads to biases of less than 0 . 05%.  \nWhile it has long been known that reconstruction reduces the nonlinear damping of the BAO in galaxy clustering due to relatively long-wavelength modes, there has been significant confusion in the literature about the appropriate way to model their effect analytically, with many works suggesting a complicated scenario where the displaced galaxies and randoms have different damping scales (e.g. Ding et al. 2018; Chen et al. 2019b). Here we show that in the RecSym scheme, where both galaxies and randoms are moved by the same redshift-space displacement, this is in fact not the case; the long-modes in fact produce a simple and unified exponential damping across both the shifted randoms and galaxies. On the other hand, we show that the RecIso scheme adopted by previous galaxy surveys like BOSS and eBOSS results in a complicated damping form where the effect of modes longer than the BAO are not controlled, leading to non-negligible shifts in the measured BAO in the process. For this reason, the RecSym scheme is the preferred one, adopted for DESI Y1 and which we recommend for future BAO analyses, though we note that in practice the differences are relatively small, especially  \nTable 1. Different contributions to the DESI BAO theory and modelling systematic error budget for 𝛼 iso and 𝛼 ap considered in this work. The section column denotes where the discussion of this contribution primarily occurs. In the theory cases, where the exact contributions depend on the nature of the tracers (galaxy bias, redshift, etc,), we have estimated a conservative value considering all DESI tracers/redshift bins.  \nat the statistical level of DESI year one and existing surveys, and we refer the reader to the accompanying refs. Chen (2024); Paillas et al. (2024) for further validation on simulations. In addition, in the process of preparing this paper, we learned that Sugiyama (2024) had reached a similar conclusion on the cancellation of long-wavelength modes in reconstruction schemes where galaxies and randoms are moved by the same reconstructed displacement. Where our paper focused on deriving the IR-resummation of the reconstructed BAO and the cancellation of the mode-coupling induced phase shift, Sugiyama (2024) instead performed a more detailed and general look at the cancellation of long-wavelength displacements in the power spectrum as a whole. We refer readers more interested in the latter aspect to that work, which we have coordinated to appear at a similar time with this one and which should be rather complementary.  \nIn the second part of the paper, we investigate the impact of a number of other modelling choices, using numerical fits to a set of highly precise simulations. We arrive at a simple set of best-practices for how to model the BAO, enumerated in §7.8. We demonstrate that if these are followed we expect a conservative systematic modelling error budget of 0 . 05% and 0 . 12% on 𝛼 iso and 𝛼 ap respectively. Any larger deviations from the expected 𝛼 values can be attributed to measurement noise. As will be demonstrated in upcoming publications, analysis of the DESI Y1 data follows these best-practices.  \nTable 1 summarises our systematic error contributions from both theory and modelling. Combining our two prongs of investigation we find that the combined systematics in BAO fitting can be controlled, conservatively, to within 0 . 1%for 𝛼 iso and 0 . 2%for 𝛼 ap . We arrive at these numbers by summing in quadrature all the detected systematics in Table 1 and rounding up to account for the unresolved systematics for which we have only upper limits. These are the values that will be adopted for modelling and theory systematics in the DESI Y1 analysis and are substantially less than the aggregate 0 . 0 < 𝑧 < 1 . 1 precision on the BAO expected from this data. They are also more than a factor of 2 less than the expected aggregate BAO precision for the full DESI sample, indicating that even with this exceptional dataset, we expect to remain in a regime where the uncertainties are dominated by statistics rather than systematics.  \nLooking to the future, it is clear that in order to further validate modeling systematics for the BAO beyond that which has been shown here, a method to test our models beyond running large numbers of N-body simulations will have to be found. Firstly, the statistical errors even from fitting sample-variance reduced DESI mocks in this work prohibit us from clearly identifying (i.e., at ≥ 3 𝜎 ) systematic  \nerrors at levels below ∼ 0 . 1%. Producing more simulations with the required precision and fidelity to overcome this will be rather costly. Secondly, we are in a realm of accuracy where the potential errors induced by galaxy clustering physics that cannot be easily introduced into N-body simulations are now important and, in addition, numerical simulations struggle due to approximations in the initial condition generation, finite computational resources and numerical approximations made during the evolution and post-processing of the simulations themselves (Angulo & Hahn 2022; Grove et al. 2022; Ding et al. 2022). While sufficient for DESI, addressing this is of the utmost importance if we aim to determine the limiting systematic precision with post-DESI BAO measurements.", 'ACKNOWLEDGEMENTS': "We thank Naonori Sugiyama for useful discussions, particularly on the perturbative treatment of infrared effects in reconstruction. We thank Alejandro Aviles and Florian Beutler for graciously agreeing to serve on the internal DESI review committee and providing helpful feedback and additionally thank Will Percival and Sesh Nadathur for useful comments.  \nSC acknowledges the support of the National Science Foundation at the Institute for Advanced Study. CH acknowledges support from the Australian Government through the Australian Research Council's Laureate Fellowship (project FL180100168) and Discovery Project (project DP20220101395) funding schemes. H-JS acknowledges support from the U.S. Department of Energy, Office of Science, Office of High Energy Physics under grant No. DESC0019091 and No. DE-SC0023241. H-JS also acknowledges support from Lawrence Berkeley National Laboratory and the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123 during the sabbatical visit. The author list was produced using the publicly available mkauthlist . 11  \nThis material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF's National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/ collaborating-institutions . Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies.  \nThe authors are honored to be permitted to conduct scientific research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.", 'DATA AVAILABILITY': "Data from the plots in this paper are available on Zenodo ( https: //zenodo.org/records/10685759 ) as part of DESI's Data Management Plan.", 'REFERENCES': "Lange J. U., 2023, MNRAS, 525, 3181  \nLazeyras T., Schmidt F., 2018, J. Cosmology Astropart. Phys., 2018, 008  \nLemaître G., 1931, MNRAS, 91, 483  \nLewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473  \nMaksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS, 508, 4017  \nMatsubara T., 2008, Phys. Rev. D, 78, 083519  \nMeiksin A., White M., Peacock J. A., 1999, MNRAS, 304, 851  \nMena-Fernández J., et al., 2024, arXiv e-prints, p. arXiv:2404.03008  \nMirbabayi M., Schmidt F., Zaldarriaga M., 2015, J. Cosmology Astropart. Phys., 2015, 030  \nMoon J., et al., 2023, MNRAS, 525, 5406  \nNoda E., Peloso M., Pietroni M., 2020, Physics of the Dark Universe, 29, 100579  \n- Noh Y., White M., Padmanabhan N., 2009, Phys. Rev. D, 80, 123501  \nPadmanabhan N., White M., 2008, Phys. Rev. D, 77, 123540  \nPadmanabhan N., White M., 2009, Phys. Rev. D, 80, 063508  \nPadmanabhan N., White M., Cohn J. D., 2009, Phys. Rev. D, 79, 063523  \n- Padmanabhan N., Xu X., Eisenstein D. J., Scalzo R., Cuesta A. J., Mehta K. T., Kazin E., 2012, MNRAS, 427, 2132  \nPaillas E., et al., 2024, arXiv e-prints, p. arXiv:2404.03005  \nPark C., Vogeley M. S., Geller M. J., Huchra J. P., 1994, ApJ, 431, 569  \nPeacock J. A., 1992, MNRAS, 258, 581  \nPeacock J. A., Dodds S. J., 1994, MNRAS, 267, 1020  \nPeebles P. J. E., 1976, Ap&SS, 45, 3  \nPeebles P. J. E., Yu J. T., 1970, ApJ, 162, 815  \nPercival W. J., et al., 2014, MNRAS, 439, 2531  \nPercival W. J., Friedrich O., Sellentin E., Heavens A., 2022, MNRAS, 510, 3207  \nPérez-Fernández A., et al., 2024, arXiv e-prints, p. arXiv:2406.06085  \nPhilcox O. H. E., Ivanov M. M., Simonović M., Zaldarriaga M., 2020, J. Cosmology Astropart. Phys., 2020, 032  \nPlanck Collaboration et al., 2020, A&A, 641, A6  \nRaccanelli A., Vlah Z., 2023, Phys. Rev. D, 108, 043537  \nRashkovetskyi M., et al., 2024, arXiv e-prints, p. arXiv:2404.03007  \nReimberg P., Bernardeau F., Pitrou C., 2016, J. Cosmology Astropart. Phys., 2016, 048  \nRobertson H. P., 1935, ApJ, 82, 284 Ross A. J., et al., 2017, MNRAS, 464, 1168 Schmidt F., 2016, Phys. Rev. D, 94, 063508  \nSchmittfull M., Feng Y., Beutler F., Sherwin B., Chu M. Y., 2015, Phys. Rev. D, 92, 123522  \nSeljak U., McDonald P., 2011, J. Cosmology Astropart. Phys., 2011, 039  \nSenatore L., Trevisan G., 2018, J. Cosmology Astropart. Phys., 2018, 019  \nSenatore L., Zaldarriaga M., 2015, J. Cosmology Astropart. Phys., 2015, 013  \nSeo H.-J., Siegel E. R., Eisenstein D. J., White M., 2008, ApJ, 686, 13  \nSeo H.-J., Beutler F., Ross A. J., Saito S., 2016, MNRAS, 460, 2453  \nSherwin B. D., Zaldarriaga M., 2012, Phys. Rev. D, 85, 103523  \nSheth R. K., 1996, MNRAS, 279, 1310  \n- Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1  \nShirasaki M., Sugiyama N. S., Takahashi R., Kitaura F.-S., 2021, Phys. Rev. D, 103, 023506  \n- Slepian Z., et al., 2018, Monthly Notices of the Royal Astronomical Society, 474, 2109\n- Smith R. E., Scoccimarro R., Sheth R. K., 2008, Phys. Rev. D, 77, 043525  \nSpeagle J. S., 2020, MNRAS, 493, 3132  \n- Sugiyama N., 2024, arXiv e-prints, p. arXiv:2402.06142\n- Sugiyama N. S., Saito S., Beutler F., Seo H.-J., 2021, MNRAS, 501, 2862  \nSunyaev R. A., Zeldovich Y. B., 1970, Ap&SS, 7, 3  \nTaylor A. N., Kitching T. D., 2010, MNRAS, 408, 865  \nTseliakhovich D., Hirata C., 2010, Phys. Rev. D, 82, 083520  \nVargas-Magaña M., et al., 2014, MNRAS, 445, 2  \nVasudevan A., Ivanov M. M., Sibiryakov S., Lesgourgues J., 2019, J. Cosmology Astropart. Phys., 2019, 037  \nVlah Z., White M., 2019, J. Cosmology Astropart. Phys., 2019, 007  \nVlah Z., Seljak U., Yat Chu M., Feng Y., 2016, J. Cosmology Astropart. Phys.,  \nWalker A. G., 1937, Proceedings of the London Mathematical Society, 42, 90  \nWallisch B., 2018, PhD thesis, University of Cambridge, UK  \nWhite M., 2001, MNRAS, 321, 1 White M., 2010, arXiv e-prints, p. arXiv:1004.0250 White M., 2015, MNRAS, 450, 3822 Xu X., Padmanabhan N., Eisenstein D. J., Mehta K. T., Cuesta A. J., 2012, MNRAS, 427, 2146 Yoo J., Seljak U., 2013, Phys. Rev. D, 88, 103520 Yoo J., Dalal N., Seljak U., 2011, J. Cosmology Astropart. Phys., 2011, 018 Zunckel C., Okouma P., Muya Kasanda S., Moodley K., Bassett B. A., 2011, Physics Letters B, 696, 433 d'Amico G., Gleyzes J., Kokron N., Markovic K., Senatore L., Zhang P., Beutler F., Gil-Marín H., 2020, J. Cosmology Astropart. Phys., 2020, 005", 'APPENDIX A: REDSHIFT EVOLUTION': 'The measured power spectrum (or correlation function) is a weighted average over the entire redshift bin, given by  \nˆ 𝑃 ( k ) = ∑︁ 𝑖 𝑤 𝑖 𝑃 ( k ( 𝑧 𝑖 ) , 𝑧 𝑖 ) . (A1)  \nThe corresponding wavenumber k varies with redshift since, if the fiducial and true cosmologies do not exactly match, the AP parameters translating between distances therein will evolve with time (see Eq. 6). Within each redshift bin, the measured power spectrum is primarily sensitive to clustering at the effective redshift 𝑧 eff = " 𝑖 𝑤 𝑖 𝑧 𝑖 at which linear trends in the time evolution cancel, i.e.  \nˆ 𝑃 ( k ) = 𝑃 ( k ( 𝑧 eff )) + 1 2 ⟨( 𝑧 -𝑧 eff ) 2 ⟩ GLYPH<16> 𝑑 2 𝑃 𝑑𝑧 2 GLYPH<17> 𝑧 = 𝑧 eff + ... (A2)  \nwhere angular brackets denote the weighted mean so that the leading correction is proportional to the variance of the redshift distribution 𝜎 2 𝑧 . This motivates the interpretation of the measured AP parameters from BAO as distances at 𝑧 eff with a theoretical error set by the width of the redshift bin.  \nLet us try to gain some heuristic understanding of the kind of errors induced in the AP parameters by the finiteness of 𝜎 𝑧 . All else being equal we can write 𝑃 = 𝑃 ( 𝛼 ∥ , 𝛼 ⊥ ) such that for example  \n𝑑𝑃 𝑑𝑧 = 𝜕𝑃 𝜕𝛼 ∥ 𝑑𝛼 ∥ 𝑑𝑧 + 𝜕𝑃 𝜕𝛼 ⊥ 𝑑𝛼 ⊥ 𝑑𝑧 .  \nThe derivatives 𝜕𝑃 / 𝜕𝛼 ∥ , ⊥ are precisely the templates for the AP parameters and, if the incorrect effective redshift were used and the linear redshift evolution not entirely cancelled, the measured distances would be off by Δ 𝛼 ∥ , ⊥ = ( 𝑑𝛼 ∥ , ⊥ / 𝑑𝑧 ) Δ 𝑧 . This would clearly be undesirable since Δ 𝑧 ∼ 0 . 1 and, within even the baseline Λ CDM model we might expect 𝑑𝛼 ∥ , ⊥ / 𝑑𝑧 ∼ ΔΩ 𝑚 / Ω 𝑚 ∼ 0 . 1, leading to percent-level systematic errors.  \nWe can repeat the same analysis for the case where 𝑧 eff were correctly determined, such that the leading error is given by the second derivative. Let us further assume that, while the first derivatives lead to phase shifts in the BAO wiggles, the second derivatives lead to only in-phase corrections that do not bias distance measurements. 12 Within this approximation we have that  \n𝑑 2 𝑃 𝑑𝑧 2 ∼ 𝜕𝑃 𝜕𝛼 ∥ 𝑑 2 𝛼 ∥ 𝑑𝑧 2 + 𝜕𝑃 𝜕𝛼 ⊥ 𝑑 2 𝛼 ⊥ 𝑑𝑧 2  \n- 12 Roughly speaking we can imagine 𝑃 𝑤 ∼ cos ( 𝛼𝑘𝑟 𝑑 ) such that the second derivative is again a cosine, now multiplied by 𝑘 2 .  \nwhich we can again directly translate into an error on the AP parameters as  \nΔ 𝛼 ∥ , ⊥ = 1 2 GLYPH<16> 𝑑 2 𝛼 ∥ , ⊥ 𝑑𝑧 2 GLYPH<17> 𝑧 = 𝑧 eff 𝜎 2 𝑧 (A3)  \nFor 𝜎 𝑧 ∼ 0 . 1 and for cosmologies wherein distances are within 10% of the fiducial one this translates to errors of less than 0 . 1%, but even in the case of larger redshift bins or models with more freedom, this shift can be exactly computed, since 𝜎 𝑧 depends only on the redshift distribution and the second derivative can be computed given an cosmology, and included when performing cosmological inference.', 'APPENDIX B: WIDE-ANGLE EFFECTS': "In the main body of the text we have operated almost exclusively within the plane-parallel approximation where the line-of-sight (LOS) can be treated as a constant vector ˆ 𝑛 . This approximation is an excellent one for most of the scales on which the BAO is measured, with corrections scaling as negative powers 𝑥 = 𝑘 𝜒 where 𝜒 is the comoving distance (Reimberg et al. 2016; Castorina & White 2018). These corrections to the plane parallel limit are due to (a) the physical LOS along which galaxies are redshifted being not parallel but radial (b) the radial evolution of the galaxy selection function 𝑑 ( 𝜒 2 ¯ 𝑛 )/ 𝑑 ln 𝜒 and (c) anisotropies incurred by the numerical implementation of the power spectrum estimator. The leading effect due to (c) is a practical rather than physical one and can be absorbed into the window function (Castorina & White 2018; Beutler et al. 2019b) - since we adopt this formalism in DESI we will focus on the first two effects (DESI Collaboration 2024f), whose leading order contributions are proportional to 𝑥 -2 .  \nSince these wide-angle effects are on large scales we can evaluate them assuming the underlying dynamics are linear, in which case the first effect (a) gives (Castorina & White 2018; Beutler et al. 2019b)  \n𝜉 0 ( 𝑠, 𝜒 ) ⊃ -GLYPH<18> 4 𝑓 2 45 𝜉 0 0 + 𝑓 ( 9 𝑏 1 + 𝑓 ) 45 𝜉 2 0 GLYPH<19> GLYPH<18> 𝑠 𝜒 GLYPH<19> 2 𝜉 2 ( 𝑠, 𝜒 ) ⊃ GLYPH<18> 4 𝑓 2 45 𝜉 0 0 + 𝑓 ( 189 𝑏 1 + 53 𝑓 ) 441 𝜉 2 0 -4 𝑓 2 245 𝜉 4 0 GLYPH<19> GLYPH<18> 𝑠 𝜒 GLYPH<19> 2  \nwhile the selection function effect gives in addition  \n𝜉 0 ( 𝑠, 𝜒 ) ⊃ 4 𝑓 2 3 𝜒 2 𝜉 0 -2 + 2 𝑓 ( 𝑏 1 -𝑓 ) 3 𝜒 𝜉 1 -1 GLYPH<18> 𝑠 𝜒 GLYPH<19> 2 𝜉 2 ( 𝑠, 𝜒 ) ⊃ -8 𝑓 2 3 𝜒 2 𝜉 2 -2 -8 𝑓 ( 5 𝑏 1 + 𝑓 ) 15 𝜒 𝜉 1 -1 GLYPH<18> 𝑠 𝜒 GLYPH<19> + 4 𝑓 2 5 𝜒 𝜉 3 -1 GLYPH<18> 𝑠 𝜒 GLYPH<19> (B1)  \nin the limit of constant number density ¯ 𝑛 . 13 In the above we have defined the generalized correlation functions  \n𝜉 ℓ 𝑛 = ∫ 𝑑𝑘 𝑘 2 2 𝜋 2 𝑘 𝑛 𝑃 ( 𝑘 ) 𝑗 ℓ ( 𝑘𝑠 ) . (B2)  \nThe contributions to the power spectrum multipoles are related to the correlation function expressions above via the usual Hankel transforms. For a realistic survey, these effects have to be averaged over line-of-sight vectors fi 𝜒 , though for the exploratory purposes of this section, we will assume a fixed 𝜒 = 𝜒 ( 𝑧 ) .  \n13 Wetake this limit in order to provide a rough estimate the size of the effect on BAO measurements, but note that for a spectroscopic survey like DESI the redshift dependence of ¯ 𝑛 can be measured to very high precision so the actual contributions could be calculated very precisely if desired.  \nIn order to estimate the bias these additional contributions cause in standard BAO fits we perform a simple Fisher forecast as we did for nonlinear effects in §3. These effects are largest at low redshift, where 𝑘 𝜒 is smallest, and when the galaxy bias is low. For this reason, let us consider the case of a sample-variance limited survey with 𝑏 1 = 1 at 𝑧 = 0 . 1, i.e. the lowest redshift boundary for the BGS sample. Even in this case, the shifts to both 𝛼 ∥ , ⊥ are uniformly less than 0 . 05% using our fiducial setup with 𝑘 min = 0 . 02 ℎ Mpc -1 . The biggest shifts are due to the selection function contribution since the odd Bessel functions are out of phase with the even ones. Nonetheless, while the wide-angle effects can reach sizes of up to several percent of the total signal at 𝑘 min their sizes are not amplified by the large BAO scale as is the nonlinear phase shift derived in §3, and furthermore since they roughly decay with wavenumber as 𝑘 -2 their overall shape is very different from a shift in the BAO scale ∝ 𝑘𝑃 ' 𝑤 ( 𝑘 ) .  \nLet us comment briefly on some other large-scale effects on the observed galaxy 2-point function not covered in the main text, such as unequal time and relativistic effects (see e.g. Grimm et al. 2020; Castorina & Di Dio 2022; Raccanelli & Vlah 2023 for recent examples). In general, these effects are suppressed on subhorizon scales like the wide-angle effects above and, moreover, do not have the exact shape of shifts in the BAO scale. For example, Raccanelli & Vlah (2023) recently showed that unequal time effects can contribute a term proportional to 𝑘𝑃 ' ( 𝑘 ) ; however, this term is suppressed on large scales and is proportional to an additional factor of ( 𝐻 ( 𝑧 )/ 𝑘 ) 2 thus it contributes negligibly to the measured BAO scale. As we do not expect any of the large-scale effects on the galaxy power spectrum to significantly bias DESI BAO measurements, we leave their detailed analysis to future work.", 'APPENDIX C: THE NONLINEAR BAO SHIFT IN REDSHIFT SPACE': "In order to compute the nonlinear BAO shift in redshift space in Eulerian perturbation theory we need to compute the galaxy density to second-order  \n𝛿 ( 2 ) 𝑔,𝑠 = ∫ p 𝑍 2 ( p , k -p ) 𝛿 lin ( p ) 𝛿 lin ( k -p )  \nwhere the second-order redshift space kernel is given in Fourier space as (see e.g. Bernardeau et al. 2002)  \n𝑍 2 ( p 1 , p 2 ) = 1 2 𝑏 2 + 𝑏 𝑠 ( p 1 · p 2 ) 2 𝑝 2 1 𝑝 2 2 -1 3 ! + 𝑏 1 𝐹 2 ( p 1 , p 2 ) + 𝑓 𝜇 2 𝐺 2 ( p 1 , p 2 ) + 𝑓 𝑘 𝜇 2 ˆ 𝑛 · p 1 𝑝 2 1 ( 𝑏 1 + 𝑓 𝜇 2 2 ) + ˆ 𝑛 · p 2 𝑝 2 2 ( 𝑏 1 + 𝑓 𝜇 2 1 ) ! (C1)  \nand we have defined 𝜇 𝑛 = ˆ 𝑛 · ˆ 𝑝 𝑛 and ˆ 𝑘 = p 1 + p 2 , with 𝜇 = ˆ 𝑛 · ˆ 𝑘 . We have also used the standard definitions of the matter density and velocity divergence kernels for 𝐹 2 and 𝐺 2 and will return to the cubic kernels 𝑍 3 when discussing BAO damping below.  \nThe largest out-of-phase contributions to the BAO in the galaxy power spectrum are sourced by the square of the quadratic contribution above (Crocce & Scoccimarro 2008; Padmanabhan & White 2009; Sherwin & Zaldarriaga 2012)  \n𝑃 ( 22 ) 𝑔,𝑠 ( k ) = 2 ∫ p 𝑍 2 ( p , k -p ) 2 𝑃 lin ( p ) 𝑃 lin ( k -p ) .  \nSpecifically, this term contains contributions from modes longer than  \nthe BAO proportional to a shift, i.e. logarithmic derivative, of the BAO template, e.g.  \n∫ | p | ≲ 1 / 𝑟 𝑑 GLYPH<18> k · p 𝑝 2 GLYPH<19> 𝑃 lin ( 𝑝 ) GLYPH<18> -k · p 𝑘 GLYPH<19> 𝑃 ' 𝑤 ( 𝑘 ) (C2)  \nwhere the first factor comes from the displacement Ψ · ∇ 𝛿 and the second one comes from expanding 𝑃 (| k -𝒒 |) for small 𝒒 . The leading terms will be those that scale as 𝑝 0 in the integrand like the one above. The bound on p comes the fact that 𝑃 (| k -𝒒 |) begins to oscillate for larger 𝒒 . Note that there are in fact two long-wavelength regions in the 22 integral since both p and k -p can be small. Additionally, the inclusion of redshift-space distortions can introduce powers of the line-of-sight angle 𝜇 into the above integral since displacements are amplified along ˆ 𝑛 .  \nIn order to compute all such terms to one-loop order in the redshiftspace galaxy power spectrum we use M/a.pc/t.pc/h.pc/e.pc/m.pc/a.pc/t.pc/i.pc/c.pc/a.pc, setting up our calculation in the spherical coordinate system  \nˆ 𝑘 = { 0 , 0 , 1 } , ˆ 𝑛 = GLYPH<26> √︃ 1 -𝜇 2 , 0 , 𝜇 GLYPH<27> ˆ 𝑝 = n√︁ 1 -𝜈 2 cos 𝜙, √︁ 1 -𝜈 2 sin 𝜙, 𝜈 o . (C3)  \nsuch that the 𝜙 and 𝜈 dependence can be performed analytically in the small p limit. Specifically, we Taylor-series expand 𝑍 2 2 and 𝑃 lin ( k -p ) in 𝑝 = | p | , perform the angular integrals in 𝜈, 𝜙 , and compute the derivative of this integral with respect to 𝑃 ' lin ( 𝑘 ) , leading to Equation 16.", 'APPENDIX D: MODELING RECONSTRUCTION IN EULERIAN PERTURBATION THEORY': 'In Appendix C we outlined the calculation of the nonlinear BAO shift in the galaxy power spectrum. This shift is, however, known to be significantly reduced by reconstruction, which we will show analytically using perturbation theory in this appendix.', 'D1 Post-Reconstruction BAO Shift in Real Space': "We begin by considering the case of matter in real space. The smoothed matter density yields displacements Ψ ( 𝑛 ) 𝑟 ( 𝒙 ) = ∇ -1 (S ∗ 𝛿 ( 𝑛 ) 𝑚 ( 𝒙 )) which can be used to compute post-reconstruction densities as  \n( 2 𝜋 ) 3 𝛿 𝐷 ( k ) + 𝛿 𝑎 ( k ) = ∫ 𝑑 3 𝒙 𝑒 -𝑖 k · 𝒙 -𝑖 k · Ψ 𝑟 ( 𝒙 ) ( 1 + 𝛿 𝑎 ( 𝒙 )) = ∞ ∑︁ 𝑛 = 0 (-𝑖 ) 𝑛 𝑛 ! FT GLYPH<2> ( k · Ψ 𝑟 ) 𝑛 ( 1 + 𝛿 𝑎 ) GLYPH<3> ( k ) . (D1)  \nHere 𝛿 𝑎 ( 𝒙 ) is the pre-reconstruction density, e.g. 𝛿 𝑑 ( 𝒙 ) = 𝛿 𝑚 ( 𝒙 ) and 𝛿 𝑠 ( 𝒙 ) = 0. We then have at second order (Hikage et al. 2017)  \n𝐾 ( 2 ) 𝑑 ( p 1 , p 2 ) = ( 1 - S k ) 𝐹 2 ( p 1 , p 2 ) -1 2 k · p 1 𝑝 2 1 S p 1 + k · p 2 𝑝 2 2 S p 2 ! + 1 2 k · p 1 𝑝 2 1 ! k · p 2 𝑝 2 2 ! S p 1 S p 2  \n𝐾 ( 2 ) 𝑠 ( p 1 , p 2 ) = -S k 𝐹 2 ( p 1 , p 2 ) + 1 2 k · p 1 𝑝 2 1 ! k · p 2 𝑝 2 2 ! S p 1 S p 2 (D2)  \nwhere we have used the shorthand S p = S( p ) .  \nWe can now proceed to compute the nonlinear BAO shift as in the pre-reconstruction case with these kernels by considering the products 𝐾 ( 2 ) 𝑎 ( p , k -p ) 𝐾 ( 2 ) 𝑏 ( p , k -p ) . If the smoothing scale of S is significantly smaller than the BAO scale 𝑟 𝑑 then for the longwavelength modes p contributing to the shift we can use S p ≈ 1 and S k -p ≈ S k . For example, performing this procedure for the autospectrum of the shifted randoms yields  \n𝑃 𝑠𝑠 ⊃ 𝑘𝑃 ' 𝑤 ( 𝑘 ) ∫ p GLYPH<18> -1 3 S 2 p S 2 k -p + 82 105 S p S k -p S k -47 105 S 2 k GLYPH<19> 𝑃 lin ( p ) ≈ 𝑘𝑃 ' 𝑤 ( 𝑘 ) ∫ p GLYPH<18> -1 3 S 2 k + 82 105 S 2 k -47 105 S 2 k GLYPH<19> 𝑃 lin ( p ) = 0 , (D3)  \ni.e. a vanishing contribution to the BAO phase shift. Computing this contribution for the 𝑑𝑑 and 𝑑𝑠 spectra similarly yields that they vanish.  \nWe can also repeat the above calculations for galaxies rather than matter. In this case we have instead Ψ ( 𝑛 ) 𝑟 ( 𝒙 ) = ∇ -1 (S∗ 𝛿 ( 𝑛 ) 𝑔 ( 𝒙 )/ 𝑏 1 ) and 𝛿 𝑑 ( 𝒙 ) = 𝛿 𝑔 ( 𝒙 ) , which we can nonetheless plug back into Equation D1 to derive 𝐾 ( 2 ) 𝑑,𝑠 for post-reconstruction galaxies. In this case, an equivalent calculation to the above yields again that the nonlinear shift is exactly cancelled by the contributions introduced by reconstruction. This is true regardless of the details of the nonlinear biasing since the nonlinear shift is a product of the bias and the quadratic shift term Ψ · ∇ 𝛿 ; the latter is always removed by reconstruction given that the correct 𝑏 1 is used. Indeed, the above calculation can be redone when  \nΨ 𝑟 = ∇ -1 S ∗ 𝛿 𝑔 ( 𝒙 ) 𝑏 fid 1 ! (D4)  \nallowing for 𝑏 fid 1 ≠ 𝑏 1 . Doing so yields a shift  \nΔ 𝛼 ∝ 𝑏 fid 1 𝑏 1 -1 ! 𝜎 2 𝑠 (D5)  \nwhere the proportionality constant is of the same order as the usual nonlinear shift pre-reconstruction, i.e. the residual shift in the BAO post-reconstruction is proportional to how well the 'true' matter overdensity can be estimated to produce the displacements. It is worth noting that the post-reconstruction shift left over due to incorrect assumptions will generally be in addition scale dependent and may not be entirely degenerate with BAO measurements, so that our estimates here should be taken to be conservative.", 'D2 Post-Reconstruction BAO Damping in Real Space': 'We can also extend the perturbation theory analysis above to look at the damping of the BAO after reconstruction. More specifically, we can use 1-loop perturbation theory to examine what long-wavelength contributions get enhanced by the BAO such that they need to be resummed.  \nBefore we do so it is useful to recall what happens prereconstruction. In this case we have for example in the matter power spectrum that (Bernardeau et al. 2002)  \n𝑃 loop 𝑚𝑚 ( k ) = ∫ p 2 𝐹 2 ( p , k -p ) 2 𝑃 lin ( 𝑝 ) 𝑃 lin (| k -p |) + 6 𝐹 3 ( k , p , -p ) 𝑃 lin ( 𝑝 ) 𝑃 lin ( 𝑘 ) .  \nBoth 𝐹 2 2 and 𝐹 3 have poles proportional to 1 / 𝑝 2 𝑛 in the longwavelength limit where one of their arguments 𝑝 𝑛 → 0, such that  \ntheir individual contributions will be quite sensitive to these IR contributions. However, if 𝑃 lin is a smooth function we can approximate 𝑃 lin (| k -p |) ≈ 𝑃 lin ( 𝑘 ) to leading order, and with some algebra it is possible to show that the poles at p , k -p → 0 in the 𝐹 2 2 piece cancel the p → 0 pole in the 𝐹 3 piece. On the other hand, if the power spectrum contains a non-smooth component such as the BAO ( ∼ sin ( 𝑟 𝑑 𝑘 ) ) then the derivative of 𝑃 lin with respect to 𝑘 is enhanced by the oscillation frequency 𝑟 𝑑 . In this case, the infrared contribution can be written as (Baldauf et al. 2015; Senatore & Trevisan 2018)  \n𝑃 loop 𝑚𝑚 ⊃ -∫ 𝑝 ≪ 𝑘 p ( k · p ) 2 𝑝 4 𝑃 lin ( 𝑝 ) ( 𝑃 𝑤 ( 𝑘 ) -𝑃 𝑤 (| k -p |)) = -∫ 𝑑 3 𝒙 𝑒 -𝑖 k · 𝒙 𝜉 𝑤 ( 𝒙 ) ∫ p , 𝑝 ≪ 𝑘 ( k · p ) 2 ) 𝑝 4 𝑃 lin ( 𝑝 ) GLYPH<16> 1 -𝑒 𝑖 p · 𝒙 GLYPH<17> ≈ -1 2 𝑘 𝑖 𝑘 𝑗 ∫ 𝑑 3 𝒙 𝑒 -𝑖 k · 𝒙 𝜉 𝑤 ( 𝒙 ) GLYPH<2> 𝐴 𝑖 𝑗 ( 𝒙 ) GLYPH<3> 𝑥 = 𝑟 𝑑 ≈ -1 2 𝑘 2 Σ 2 NL 𝑃 𝑤 ( 𝑘 ) . (D6)  \nwhere 𝜉 𝑤 is the Fourier transform of the 𝑃 𝑤 and we have used that 𝜉 𝑤 is localized at 𝑟 𝑑 to perform a saddle point approximation between the second and third lines. This is none other than the 1-loop contribution to the BAO damping 𝑒 -1 2 𝑘 2 Σ 2 NL 𝑃 𝑤 and uses the same saddle-point approximation in LPT described in §3, with the result that the BAO is damped by (numerically) large contributions from long-wavelength modes.  \nWe can perform the same calculation post-reconstruction. For brevity, we will only consider the case of the shifted randoms in matter reconstruction in order to highlight why our result in §4 differs from the existing literature-the calculation for shifted galaxies etc. is entirely analogous. In order to do so we need to compute the third-order contribution in Equation D1  \n𝐾 ( 3 ) 𝑠 ( p 1 , p 2 , p 3 ) = -S k 𝐹 3 ( p 1 , p 2 , p 3 ) + 1 3 ( k · p 1 )( k · p 23 ) 𝑝 2 1 𝑝 2 23 S 1 S 23 + et . cycl . ! -1 6 ( k · p 1 )( k · p 2 )( k · p 3 ) 𝑝 2 1 𝑝 2 2 𝑝 2 3 ! S 1 S 2 S 3  \nwhere we have used numerical indices to indicate arguments by p 𝑛 and used the shorthand p 𝑛𝑚 = p 𝑛 + p 𝑚 . Plugging this into the 1-loop power spectrum then yields  \n𝑃  \nloop 𝑠𝑠 = ∫ p 2 𝐾 ( 2 ) 𝑠 ( p , k -p ) 2 𝑃 lin ( 𝑝 ) 𝑃 lin (| k -p |) + 6 𝐾 ( 3 ) 𝑠 ( k , p , -p ) 𝑃 lin ( 𝑝 ) 𝑃 lin ( 𝑘 ) ⊃ -S 2 k ∫ 𝑝 ≪ 𝑘 p ( k · p ) 2 𝑝 4 ( 1 - S p ) 2 𝑃 lin ( 𝑝 ) ( 𝑃 𝑤 ( 𝑘 ) -𝑃 𝑤 (| k -p |)) . (D7)  \nThis is identical to the integral except that (1) there is a factor of S 2 k in front, since 𝑃 𝑠𝑠 is given by S 2 𝑃 lin at linear order and more (2) the integrand is suppressed by ( 1 - S p ) 2 internally. This is because, as described in §4, the shifted field is in fact displaced by bulk motions given to leading order by ( 1 - S) ∗ Ψ even though the randoms themselves are only moved by -S ∗ Ψ . This is also reflected by the fact that the nonlinear BAO shift derived above cancels independently in 𝑃 𝑠𝑠 by virtue of the Ψ · ∇ 𝛿 contribution cancelling there as well. Similar calculations show that the same result occurs in both 𝑃 𝑑𝑑,𝑑𝑠 as well, such that all of the pieces of the reconstructed power spectrum  \nhave the same BAO damping. The case for galaxies is also entirely similar.', 'D3 Nonlinear BAO in Redshift Space': "Finally, let us extend our treatment of the nonlinear BAO postreconstruction at 1-loop in EPT to include redshift-space distortions. The order 𝑛 contribution to the reconstructed displacement is given by  \nΨ ( 𝑛 ) 𝑟,𝑠 ( k ) = -GLYPH<18> 𝑖 k 𝑘 2 GLYPH<19> ' › 𝛿 ( 𝑛 ) 𝑔,𝑠 ( k ) 𝑏 1 + 𝑓 𝜇 2 ' fi . (D8)  \n«  \n‹  \nIn RecSym both the galaxies and randoms are shifted by this displacement multiplied by the redshift-space transformation matrix 𝑅 𝑖 𝑗 = 𝛿 𝑖 𝑗 + 𝑓 ˆ 𝑛 𝑖 ˆ 𝑛 𝑗 . Let us again consider the case of matter such that 𝑏 1 = 1 for simplicity. In this case, the quadratic kernel for the shifted randoms is now  \n𝑍 ( 2 ) 𝑠 ( p 1 , p 2 ) = (-𝑖𝑘 𝑖 ) 𝑅 𝑖 𝑗 -GLYPH<18> 𝑖𝑘 𝑗 𝑘 2 GLYPH<19> 𝑍 ( 2 ) ( p 1 , p 2 ) 1 + 𝑓 𝜇 2 !! S k -1 2 𝑘 𝑖 𝑅 𝑖 𝑗 -𝑖 𝑝 1 , 𝑗 𝑝 2 1 𝑍 1 ( p 1 ) 1 + 𝑓 𝜇 2 1 !! 𝑘 𝑖 𝑅 𝑖 𝑗 -𝑖 𝑝 2 , 𝑗 𝑝 2 2 𝑍 1 ( p 2 ) 1 + 𝑓 𝜇 2 2 !! S 1 S 2 (D9)  \nwhere as usual k = p 1 + p 2 and 𝜇 𝑛 = ˆ 𝑛 · ˆ 𝑝 𝑛 . We have included the dependence on 𝑍 1 explicitly but note that 𝑍 1 ( p 𝑛 )/( 1 + 𝑓 𝜇 2 𝑛 ) = 1. Repeating the real-space calculation using the coordinate system of Equation C3, and explicitly integrating over the angular dependence of the long mode p , it is straightforward, if laborious, to show that the nonlinear shift proportional to 𝜎 2 𝑠 𝑘𝑃 ' 𝑤 ( 𝑘 ) vanishes. In addition, we can also see that the arguments about the damping scale being set by ( 1 -S)∗ Ψ also hold in redshift space since the relevant (resummed) contribution comes from powers of the redshift-space displacement sourced by S p 𝑛 𝛿 ( 1 ) ( p 𝑛 )/( 𝑏 1 + 𝑓 𝜇 2 𝑛 ) . These smooth displacements cancel the exponentiated Zeldovich displacements on large scales at all orders, such that their total contribution can be resummed as claimed in the main text. The case of displaced galaxies is entirely analogous.  \nFinally, let us repeat the calculation for RecIso . In RecIso the particles are shifted by the displacement without multiplying by 𝑅 𝑖 𝑗 , and we can reproduce this in our calculations by taking 𝑅 𝑖 𝑗 → 𝛿 𝑖 𝑗 in the above. Performing the same calculation yields a BAO shift in RecIso of  \n' › › 7 𝑓 3 GLYPH<16> 4 -9 𝜇 2 GLYPH<17> 𝜇 4 + 𝑓 2 GLYPH<16> 14 -73 𝜇 2 GLYPH<17> 𝜇 2 -12 𝑓 𝜇 2 105 GLYPH<0> 1 + 𝑓 𝜇 2 GLYPH<1> ' fi fi 𝜎 2 𝑠 𝑊 2 ( 𝑘 ) 𝑘𝑃 ' 𝑤 ( 𝑘 )  \n«  \n‹  \nrepresenting a nonzero shift unlike RecSym .", 'APPENDIX E: HANKEL TRANSFORM OF BROADBAND SPLINE MODEL': 'The two relevant Hankel-transformed spline coefficients for the correlation function were evaluated to  \n𝐵 2 , 0 ( 𝑥 ) = 2 𝑥 6 GLYPH<16> 𝑥 3 Si ( 𝑥 ) -2 𝑥 3 Si ( 2 𝑥 ) + 𝑥 2 cos ( 𝑥 ) -𝑥 2 cos ( 2 𝑥 ) -𝑥 sin ( 𝑥 ) -16 cos ( 𝑥 ) + 4 cos ( 2 𝑥 ) + 𝑥 sin ( 𝑥 ) cos ( 𝑥 ) + 12 GLYPH<17>  \n𝐵  \n2  \n,  \n1  \n(  \n𝑥  \n)  \n=  \n-  \n1  \n2  \n𝑥  \n3  \nGLYPH<16>  \n6  \n𝑥  \n+  \n-  \n8  \n𝑥  \n6  \n𝑥  \nsin  \n) -  \n𝑥  \n+  \n6  \n𝑥  \n(  \n𝑥  \n32  \n𝑥  \n2  \ncos  \n) +  \n8  \n𝑥  \nsin  \nSi  \n(  \n𝑥  \n(  \n2  \n𝑥  \n) -  \n(  \n) +  \n16  \n𝑥  \n2  \n𝑥  \n) -  \n27  \n𝑥  \n2  \nSi  \ncos  \n3  \n𝑥  \nsin  \n(  \n2  \n𝑥  \n(  \n3  \n𝑥  \n)  \n) +  \n3  \n𝑥  \n)  \n-96 cos ( 𝑥 ) + 64 cos ( 2 𝑥 ) -16 cos ( 3 𝑥 ) + 48 GLYPH<17> (E1)  \nwhere Si ( 𝑥 ) is the sine integral ∫ 𝑧 0 𝑑𝑡 sin ( 𝑡 )/ 𝑡 .', 'APPENDIX F: ANALYTIC MARGINALISATION OF THE BROADBAND': 'Our full model for the BAO contains a large number of nuisance parameters designed to recover unbiased constraints on the two parameters we actually care about ( 𝛼 iso and 𝛼 ap ) in the presence of other physics that affects the shape of the galaxy clustering, survey selection effects and measurement noise.  \nIn our base model for the galaxy clustering however, most of these nuisance parameters enter at linear order, and so can be marginalised over analytically at the level of the likelihood without needing to be sampled, expensively and numerically, as part of an MCMC or Nested sampling algorithm (Bridle et al. 2002; Taylor & Kitching 2010; d\'Amico et al. 2020). This is useful as we are typically not interested in the constraints on these other parameters.  \nTo demonstrate how this works, we rewrite the multipole version of our clustering model from Eq. 5 into two parts, where the first combines the wiggle and no-wiggle model components, and the latter is the broadband, making clear its dependence on the linear nuisance parameters 𝑎 ℓ,𝑛  \n𝑃 𝑔𝑔,ℓ ( 𝑘 ) = ˜ 𝑃 𝑔𝑔,ℓ ( 𝑘 ) + ∑︁ 𝑛 𝑎 ℓ,𝑛 𝜕 D ℓ ( 𝑘 ) 𝜕𝑎 ℓ,𝑛 . (F1)  \nWhen performing a BAO fit, this model is then convolved with the survey window function (matrix W ℓℓ \' ) to give a new vector 𝑷 𝑊 ℓ = W ℓℓ \' 𝑷 𝑔𝑔,ℓ \' = ˜ 𝑷 𝑊 ℓ + " 𝑛 𝑎 ℓ,𝑛 𝑫 𝑊 ℓ,𝑛 where we have still kept the windowed model in two parts, those that do and do not depend on the linear order nuisance parameters respectively. The matrix 𝐷 𝑊 ℓ,𝑛 then gives a set of purely linear templates in the power spectrum model. These templates can in general depend on the nonlinear degrees of freedom (e.g. bias and dilation parameters) though they do not in our particular case.  \nConcatenating the different multipoles, we then take the common approach of using a Gaussian likelihood to describe our data ( ˆ 𝑷 ) given the model with parameters 𝜽 ,  \nL( ˆ 𝑷 | 𝜽 ) ∝ Exp GLYPH<20> -1 2 [ ˆ 𝑷 -𝑷 𝑊 ] C -1 [ ˆ 𝑷 -𝑷 𝑊 ] 𝑇 GLYPH<21> . (F2)  \nThe act of marginalising over the nuisance parameters, before or after sampling them, is akin to performing the integral  \nL( ˆ 𝑷 | 𝛀 ) ∝ ∫ 𝑑 fi 𝒂 L( ˆ 𝑷 | 𝜽 ) 𝑃 (fi 𝒂 ) (F3)  \nwhere 𝑃 (fi 𝒂 ) are the prior probabilities for the various 𝑎 ℓ,𝑛 and 𝛀 denotes all the other free parameters of the model. When using the Gaussian likelihood in Eq. F2 and a flat (unbounded) prior for 𝑃 (fi 𝒂 ) , this integral has an exact analytic solution,  \nL( ˆ 𝑷 | 𝛀 ) ∝ 1 | F | 1 / 2 Exp GLYPH<20> -1 2 [ ˆ 𝑷 -˜ 𝑷 𝑊 ] ˜ C -1 [ ˆ 𝑷 -˜ 𝑷 𝑊 ] 𝑇 GLYPH<21> , (F4)  \n(  \n3  \nSi  \n(  \n2  \n3  \n3  \n9  \n𝑥  \n2  \nwhere  \n(  \n)  \nF = D 𝑊 C -1 D 𝑊,𝑇 (F5)  \nis a square matrix with length equal to the number of nuisance parameters, and  \n˜ C -1 = C -1 -C -1 D 𝑊,𝑇 F -1 D 𝑊 C -1 . (F6)  \nThis has the interpretation that the full covariance is given by the data covariance plus a theory covariance, such that the total precision matrix is reduced. Similar analytic solutions can be obtained for the cases where 𝑃 (fi 𝒂 ) takes a Gaussian, in which cases one simply replaces  \nF → F + Σ (F7)  \nwhere Σ is the covariance matrix for the Gaussian prior on fi 𝒂 .  \nGiven the above expressions, one can also identify the best-fit values of the linear order nuisance parameters given the remainder of the model ˜ 𝑷 𝑊 ,  \n𝒂 best -fit = F -1 D 𝑊 ˜ C -1 [ ˆ 𝑷 -˜ 𝑷 𝑊 ] 𝑇 . (F8)  \nThe improvement in 𝜒 2 obtained by using these best fit parameters is equal to the reduction in the precision matrix above.  \nIf the linear templates in D 𝑊 are independent of the nonlinear degrees of freedom, as is the case for the BAO broadband, an additional simplification occurs. In this case F and ˜ C -1 can be precomputed and so evaluation of the likelihood takes no additional time compared to the unmarginalised case. Combining with the above discussion this allows one to also take the \'partially\' marginalised approach - evaluating the best-fit broadband parameters at each sampling iteration, multiplying by D 𝑊 and adding them to the rest of the model to evaluate the standard 𝜒 2 . This is the approach taken in some past literature (Ross et al. 2017; Bautista et al. 2021) and gives the same marginalised posterior on the parameters of interest in this special case. However, as can be seen from the description above, it requires more numerical operations per likelihood evaluation and so we argue the direct marginalisation is preferable.  \nIn either case, the benefits to sampling and run-time are substantial when including analytic marginalisation. For our default model, when fitting the LRG mocks the number of free parameters in the sampling is reduced from 21 to 7. As such, we adopt full analytic marginalisation over the broadband parameters as our default.  \nAnalytic marginalisation can be carried out similarly for the correlation function, and further inspection shows that the galaxy bias (squared) could also be analytically marginalised over if desired. However, we do not include this latter option in our default methodology as tests on noisier data suggest that allowing the squared galaxy bias to go negative (as must be allowed in the case of a flat, unbounded prior) can give incorrect fits to the BAO in the correlation function, where the \'dip\' in the correlation function on scales smaller than the BAOismistaken for (a flipped version of) the BAO itself. The method of analytic marginalisation does not lend itself to mixed priors, and so an alternative for further investigation would be to investigate Gaussian priors on the galaxy bias and broadband terms.', 'APPENDIX G: AFFILIATIONS': "- 1 Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA  \ncos  \n3  \n𝑥  \n- 3 Department of Physics, University of California, Berkeley, 366 LeConte Hall MC 7300, Berkeley, CA 94720-7300, USA\n- 4 Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA\n- 5 Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA\n- Department of Physics & Astronomy, Ohio University, Athens,\n- 6 OH 45701, USA\n- 7 Physics Department, Yale University, P.O. Box 208120, New Haven, CT 06511, USA\n- 8 Physics Dept., Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA\n- 9 Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India\n- 10 University of Michigan, Ann Arbor, MI 48109, USA\n- 11 Leinweber Center for Theoretical Physics, University of Michigan, 450 Church Street, Ann Arbor, Michigan 48109-1040, USA\n- 12 NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA\n- 13 Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK\n- 14 Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK\n- 15 Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA\n- 16 Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México C.P. 04510, México\n- 17 Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China\n- 18 University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA 94720, USA\n- 19 Institut de Física dÕAltes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra Barcelona, Spain\n- 20 Ecole Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland\n- 21 Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia\n- 22 Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia\n- 23 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA\n- 24 Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA\n- 25 Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain\n- 26 Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, UK\n- 27 Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans s/n, 08913 Bellaterra, Barcelona, Spain\n- 28 Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA\n- 29 The Ohio State University, Columbus, 43210 OH, USA\n- 30 Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Spain\n- 31 Laboratoire de Physique Subatomique et de Cosmologie, 53 Avenue des Martyrs, 38000 Grenoble, France\n- 32 Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, 08010 Barcelona, Spain\n- 33 Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada\n- 34 Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada  \n- 35 IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France\n- 36 Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada\n- 37 Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 1, 85748 Garching, Germany\n- 38 Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain\n- 39 Department of Physics, Kansas State University, 116 Cardwell Hall, Manhattan, KS 66506, USA\n- 40 Department of Physics and Astronomy, Sejong University, Seoul, 143-747, Korea\n- 41 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia 42 CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain\n- 43 SLAC National Accelerator Laboratory, Menlo Park, CA 94305, USA\n- 44 National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd., Chaoyang District, Beĳing, 100012, P.R. China  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author."}