Compressed baryon acoustic oscillation analysis is robust to modified-gravity models

We study the robustness of the baryon acoustic oscillation (BAO) analysis to the underlying cosmological model. We focus on testing the standard BAO analysis that relies on the use of a template. These templates are constructed assuming a fixed fiducial cosmological model and used to extract the location of the acoustic peaks. Such “compressed analysis” had been shown to be unbiased when applied to the ΛCDM model and some of its extensions. However, it has not been known whether this type of analysis introduces biases in a wider range of cosmological models where the template may not fully capture relevant features in the BAO signal. In this study, we apply the compressed analysis to noiseless mock power spectra that are based on Horndeski models, a broad class of modified-gravity theories specified with eight additional free parameters. We study the precision and accuracy of the BAO peak-location extraction assuming DESI, DESI II, and MegaMapper survey specifications. We find that the bias in the extracted peak locations is negligible; for example, it is less than 10% of the statistical error for even the proposed future MegaMapper survey. Our findings indicate that the compressed BAO analysis is remarkably robust to the underlying cosmological model.


Introduction
Baryon acoustic oscillations (BAO) have by now established themselves as a main probe of cosmology, providing constraints on dark energy and the expansion history of the universe.The physics of the BAO has been well understood starting from the pioneering work by [1] and [2]: primordial sound waves in the baryon-photon fluid prior to recombination imprint a specific feature in the distribution of overdensities in the universe.This physical feature -the sound horizon of about r d ≈ 150 Mpc in the standard cosmological model -can be observed today in the distribution of galaxies as the scale at which there is a ∼10% excess probability for clustering.The sound horizon allows for precise measurements of the angular diameter distance and the Hubble parameter at low redshifts (z ∼ 1) where tracers of the large-scale structure -galaxies and quasars -are typically observed.The BAO feature was first detected and used to constrain cosmology nearly two decades ago [3,4].Subsequent analyses have spearheaded an increasingly effective use of the BAO feature to constrain the cosmological parameters and models (e.g., [5][6][7][8]).Because the BAO features (in Fourier space, or a single feature in configuration space) reside in the linear-clustering regime, BAO are relatively free from systematic errors associated with nonlinear physics.BAO are thus a powerful tool to constrain different cosmological models of dark matter and dark energy.
The standard BAO analysis -the one that had been most commonly applied to datafocuses on extracting the angular features corresponding to the sound horizon, while simply fitting out the broadband power spectrum.In other words, this kind of analysis fits a template that had been created assuming a fiducial cosmology.Galaxy clustering data are fitted to this template to extract the BAO feature(s), after marginalizing over many nuisance parameters which account for differences between the template and the measured broadband clustering.This standard analysis that makes use of a template is sometimes referred to as the compressed analysis (since it compresses the clustering information into the transverse and radial location of the BAO peak), and we describe it in detail in Sec. 3.This procedure has been the basis for deriving cosmology from BAO measurements starting from the earliest analyses (e.g., [3,[5][6][7][9][10][11][12]).Alternatives to the standard analysis include the so-called "direct fit" (sometimes also called "full-shape modeling"), which fully models the broadband power spectrum including the BAO peak [13][14][15][16], as well as the ShapeFit method, which is similar to the compressed analysis but includes a single additional parameter that extracts additional information about the slope of the power spectrum [17].Additionally, one can also extract information from BAO signal using the "linear point" in the correlation function instead of the BAO peaks, which has the advantage of being less sensitive to nonlinear effects [18].
It has been demonstrated that the compressed analysis is robust when one assumes the standard ΛCDM cosmological model [19,20].In other words, the compressed analysis, which extracts the BAO peak location (or rather the relative location in the transverse direction, α ⊥ and that in the radial direction, α ∥ ; we will introduce these parameters in Sec.3), recovers the true values of the cosmological parameters.This is not too surprising in the ΛCDM model, essentially because the compressed analysis is based on a template that had been constructed assuming ΛCDM.However, it is possible that models beyond ΛCDM add features to the power spectrum that cannot be well modeled by the template, and thus introduce unaccounted-for systematic errors that would bias the measured parameters.
Relatively small deviations from ΛCDM (for example wCDM model, which adds as a free parameter the dark energy equation of state w) are expected to remain robust under standard compressed BAO analysis.This has been validated to some extent by previous studies that confirmed the flexibility and effectiveness of standard BAO analysis with different methodology choices and against data generated with different cosmological models.For instance, [21] found that the extracted BAO peak-location parameters have negligible dependence on the assumed fiducial cosmologies, but their errors have a non-trivial increase when the fiducial cosmologies deviate from the test models.Similarly, [22] simulated the BAO compressed analysis assuming that the cosmological models with modified perturbations before recombination, and found no significant shifts in the extracted cosmological parameters.
However, the aforementioned studies have tested individual cosmological models, and no attempt to "sweep" through the much larger space of beyond-standard cosmological models has, to our knowledge, been attempted.The question remains therefore of precisely how robust is the compressed analysis in the presence of very general cosmological models, and how often (if at all) the standard BAO analysis fails under such models.We endeavor to answer this question, and reproduce the compressed BAO analysis algorithm, then apply it to a wide range of modified-gravity models.We will focus on Horndeski models (introduced and explained in Sec. 2) for the following two reasons: (1) this is a very general class of modifications of gravity and is arguably well motivated to potentially explain the accelerating expansion of the universe, and (2) Horndeski models have been thoroughly studied in the literature, and in particular there exist publicly available Einstein-Boltzmann codes that produce the basic cosmological observables (like the primordial matter power spectrum) for an arbitrary model from this class.
The outline of the paper is as follows.We introduce the Horndeski models and describe their parameterization that we adopt along with theory parameter priors, in Sec. 2. In Sec. 3, we describe the compressed BAO methodology, and specifically our implementation of it along with all relevant details and assumptions.In Sec. 4, we present the result of our tests for the robustness of BAO standard analysis in Horndeski models.We conclude and discuss other lines of current and future work in Sec. 5.

Horndeski models of modified gravity
We now provide a concise overview of the Horndeski models.We present the selected parameters for the test models and discuss the impact of modified gravity on the resulting matter power spectra.
[23] introduced models with the most general second-order Euler-Lagrange equations that can be obtained from the metric g µν , the scalar field ϕ, and their derivatives in fourdimensional space.Long after it was first proposed, the importance of the Horndeski framework was revisited and recognized by [24] who reduced the original Lagrangian to a combination of four base Lagrangians.In this study, we follow the Effective Field Theory (EFT) approach [25,26] that parameterizes the Horndeski models with a small number of free functions, which can be further reduced to a few parameters that control the cosmological background and perturbations.
The action in unitary gauge for the EFT of dark energy can be written as the following (e.g., [26][27][28]): where M PL is the Planck mass, δg 00 is defined as g 00 + 1, δK ν µ is the perturbation of the extrinsic curvature, δK is its trace, R is the Ricci scalar, δR (3) is the perturbation of the spatial component of the Ricci scalar, and S m (g µν , Ψ m ) is the action of matter field.There are a number of free functions here: Ω MG (t), Λ(t), c(t), M 2 (t), M1 (t), M2 (t), M3 (t), M (t), and m 2 (t).The first three functions determine the background evolution; because c(t) and Λ(t) are subject to constraints from the energy density and pressure respectively in the Friedmann equations, the background evolution in modified gravity is controlled by the single function Ω MG (t).In the literature, Ω MG (t) is often referred to as Ω(t); we use Ω MG (t) to avoid confusion with an energy-density parameter.The remaining free functions determine the evolution of perturbations.
For convenience, we redefine the second-order free functions in a dimensionless form (see [29] for an alternative parameterization).The dimensionless functions are In linearized Horndeski theory, the free functions governing the evolution of perturbations are subject to the following constraints: (2. 3) The constraints in Eq. (2.3) imply 2γ 5 = γ 3 = −γ 4 and γ 6 = 0. We adopt the following ansatz for the time-dependence of the remaining gammas: since this functional form is simple yet reasonably flexible.Similarly, we choose Ω MG (t) to have a form Ω MG (a) = Ω MG,0 (a) a s,0 . (2.5) Thus, there are eight free parameters of the Horndeski models We adopt the EFTCAMB code [30] to produce cosmological observables with Horndeski models described with the parameterization above.The Horndeski parameters above specify the perturbations, but not the background.For the latter, we adopt the ΛCDM expansion history in a flat universe, with the single free parameter Ω M = 1 − Ω Λ .
We now discuss the priors that we give to the Horndeski parameters; the priors are similar to (but not the same as) those adopted in [31].The priors on γ 1,0 and γ 2,0 are chosen based on the preferred values from current cosmological data [32].The parameter γ 3 relates the speed of gravitational waves to the speed of light via where c T is the speed of gravitational waves.Here we choose γ 3,0 = 0; as gravitational waves propagate at the speed of light.In particular, the gravitational-wave event GW170817 ruled out all Horndeski models with γ 3,0 ̸ = 0. Note that theories beyond general relativity, including the Horndeski class, allow non-luminal gravitational-wave speed at low energies.We set γ 3,0 = 0 in order to prevent models with non-luminal tensor speed at z = 0 (see [33] for further discussion).All of these priors are summarized in Table 1.The parameterization of the cosmological parameters that control the background (and their associated priors) will be discussed below, in Sec.3.1.Additionally, we impose physical stability conditions, mathematical stability conditions, and EFT additional conditions (see section IV F in [30] for details).The physical stability, including both ghost and gradient stability conditions, ensures that the background evolution is stable (see Eqs. ( 42)-( 51) in [30]).Ghost instability refers to a wrong sign of the kinetic term.Gradient instability is typically associated with a negative square of the sound speed, c 2 s < 0 in the equations of motion of perturbations, leading to unbounded growth of small-scale perturbations.Mathematical stability conditions necessitate a well-defined π-field equation, the absence of fast exponential growth in the π-field perturbations, and a well-defined equation for tensor perturbations (see Eq. ( 52) in [30]).The mathematical stability conditions ensure that the perturbation in the dark section is stable (see Eqs. ( 30)-( 32) and ( 41)-( 52) in [30] for details of the physical and mathematical stability conditions effects on the parameters Ω MG (a) and γ i (a)).The EFT additional conditions require that w(a) ≤ −1/3, which is already satisfied as we have fixed w = −1.

Simulating and measuring the BAO scale in Horndeski models
We now describe the methods that we used to test the bias of the BAO standard analysis.The procedure involves the following steps: 1. We generate mock power spectra (and their multipoles), along with their corresponding covariance, based on each assumed galaxy survey and the underlying cosmological parameters.
2. We fit the mock power spectrum multipoles using a template, thus jointly constraining about 15 cosmological and nuisance parameters.The parameters of our interest are the α parameters that describe the BAO location.
3. We quantify the bias of the test model in the standard BAO analysis utilizing the best-fitted parameters and their confidence intervals.
The first two steps in this procedure are quantified in the rest of this Section.The fourth step constitutes our principal results, as outlined in Sec. 4.

Cosmological model parameters
To scan through a range of Horndeski cosmological models, we vary the Horndeski parameters given in Table 1, as well as the cosmological parameters that specify the ΛCDM background.The Horndeski parameters p Horn are specified in Eq. (2.6), while the cosmological parameters are where h is the Hubble constant in units of 100km/s/M pc, Ω cdm h 2 and Ω b h 2 are respectively the physical cold dark matter and baryon energy densities, A s and n s are the amplitude and spectral index of primordial density fluctuations, and τ reio is the optical depth of reionization.The Horndeski parameters are sampled from flat priors listed in Table 1.In contrast, we choose a more complicated (correlated) prior for the base cosmological parameters from Eq. (3.1) in order to make them in reasonably good agreement with respect to the current data; we do so since we do not wish to study models that are obviously ruled out.We choose the base cosmological parameters that generate ΛCDM models that are within 5σ confidence interval (∆χ2 < 39.4 in six-dimensional space from Eq. (3.1)) relative to the best-fit model from Planck TTTEEE + low E data [34] 1 We illustrate this in Fig. 1, where we show the sampled ΛCDM parameters and the 5σ contours of Planck data.We fix curvature to zero, and the neutrino density to Ω ν h 2 = 0.001.
For each set of base cosmological and Horndeski parameters, we generate a matter power spectrum using in EFTCAMB.Subsequently, we utilize this derived matter power spectrum to compute P (k, µ) data following the fitting template formulation that we describe in Sec.3.3 below.We adopt values from the fiducial cosmology for the template nuisance parameters (that is, all parameters in the fit other than the relative location of the BAO peak, α ∥ and α ⊥ ; see also section 3.3).We also need to specify a fixed set of cosmological parameters that describes the BAO template.We select the parameters that are close, but not identical, to the Planck best fit: h = 0.6736, Ω cdm h 2 = 0.119, Ω b h 2 = 0.022, A s = 2.1 × 10 −9 , n s = 0.9649, and τ reio = 0.0544.
We next describe the procedure for calculating matter power spectra of Horndeski models.

Matter power spectra in Horndeski models
We use EFTCAMB [30] to generate matter power spectra predicted by Horndeski models.Note that the Horndeski parameters, listed in Eq. (2.6), only affect the perturbations and not the cosmological background (distances and volumes) 2 , instead only affecting the shape of the BAO features in Fourier space.In contrast, the base cosmological parameters do shift the BAO peak; for example, the late-universe energy densities of dark matter and dark energy control the distance to the galaxies/quasars, and hence the angular extent of the sound horizon observed at the corresponding redshift.Therefore, each one of our models has a different BAO scale than that predicted by the fiducial ΛCDM model at that redshift.While the Horndeski parameter variations by themselves do not shift the BAO, they change other features of P (k) which are degenerate with those induced by varying the base ΛCDM parameters.Thus, the overall shift of the peak location in our models is more complex than that in vanilla ΛCDM.
Fig. 2 illustrates the shifts in the location of BAO peaks for a representative sample of Horndeski models that includes variations of both the background (cosmological) and perturbation (Horndeski) parameters.Note that the amplitude of the power spectrum and the locations of BAO peaks both vary in a complex way.We illustrate the P (k) changes in more detail in Appendix A, where Fig. 5 and Fig. 6 illustrate the change of matter power spectrum when the perturbation parameters alone are varied one at a time.We also present a comparison of the best-fitted parameters for the power spectrum in Horndeski models against those in the fiducial cosmology in Table 4.We observe that the perturbation parameters alone induce some variation in the amplitude of the power spectrum, but not in the BAO location.At the same time, the base cosmological parameters do change the BAO peak locations as expected.Therefore, any given Horndeski model will have different BAO peak locations along with changes in the amplitude of the power spectrum that are potentially different than that expected from parameter variations in the ΛCDM model.This fact motivates our investigation, which is to see whether extraction of the BAO peak information in modifiedgravity models that use a fixed template that is centered around ΛCDM can recover unbiased cosmological results.

Template for the anisotropic power spectrum
We now review the standard BAO analysis -how to isolate and measure the BAO signal from our mock realizations in both isotropic and anisotropic cases.Our analysis is anisotropic, i.e., Figure 2: Twenty-four power spectra were randomly selected from our sampler, where both the Horndeski parameters and the base ΛCDM cosmological parameters were randomly sampled.The top panel shows the power spectra, the middle panel shows their ratios to the smooth fiducial power spectrum, P (k)/P smoooth, fid (k), while the lower panel shows the ratios relative to the smoothed portion of each corresponding spectrum, O lin = P (k)/P smoooth (k).
In each panel, the power spectrum evaluated in the fiducial cosmological model is plotted in red line.Note the non-negligible shifts in the BAO wiggle position.
separates the transverse and radial modes on the sky.Nevertheless, in the interest of pedagogy, we first review the isotropic analysis in order to introduce some key (and by now standard) tools.
To extract the BAO peak locations from data, it is economical to first fit a template to the power spectrum.The template assumes some fixed fiducial cosmological model and has a key feature of allowing freedom in the horizontal shift of the BAO features in k-space.Specifically (and still assuming the isotropic BAO case for the moment), the BAO shift is controlled by the α, defined as Here we have defined a generalized distance [3] where D A (z) is the angular diameter distance, H(z) is the Hubble parameter, and r d is the comoving sound horizon at the drag epoch.Here, D V is the distance that quantifies the average of the distances measured along, and perpendicular to, the line of sight to the observer.Moreover, the subscript "fid" refers to the corresponding values at the (fixed) fiducial cosmology, while D V , H(z), and r d are evaluated in the cosmological model that is being tested.
The other parameters that enter the template, to which we will refer as the "nuisance parameters", also carry potentially useful cosmological information (about e.g., the amplitude and shape of the primordial power spectrum), but are less robust than the α parameter as they are degenerate with systematic and astrophysical parameters, for example, the galaxy bias.We now introduce these remaining template parameters.We model the isotropic power spectrum following [6] where and Here B p , Σ 2 nl and A i (with i ranging from 1 to 5) are all nuisance parameters: B p accounts for potential large-scale bias, A i accounts for the possibility that P (k) sm does not match the actual data, and Σ nl characterizes the damping of BAO.Next, P (k) lin is the linear matter power spectrum, while P (k) sm,lin refers to the smooth part of the linear matter power spectrum, i.e. one without the BAO features (here we adopt a smoothing method in the configuration-space; see 3.5 for details).Now we generalize the template to allow for anisotropy in the power spectrum (and the BAO peak locations).Instead of the single parameter α, we now have α ∥ and α ⊥ that describe the BAO features in the parallel and transverse directions to the line-of-sight, respectively.They are defined as The fiducial cosmology is used to convert the measured redshift to the distance and establish the conversion factors between the fitting template and true values of r d .The values of α ∥ and α ⊥ in the fiducial cosmology are unity.We next need to link the wavenumber utilized by EFTCAMB, k ≡ |k|, and shift it to the wavenumber(s) that describe the anisotropic power spectrum in an arbitrary cosmological model.Starting with some k that we provide to the numerical code, we first consider the cosine of the angle between this wavenumber and its projection along the line of sight, µ ≡ k ∥ /k; in this way, we get the wavenumber components that are respectively parallel and perpendicular to the line of sight, k ∥ and k ⊥ .Next, we track how these two components scale to reflect the shift of the BAO peak in an arbitrary cosmological model.The "observed" wavenumbers in a given model are k aniso where k ∥ and k ⊥ are the values in the fiducial model.Finally, we can then express the coordinates of the anisotropic fitting in terms of the k and µ using the following relationships: The fitting template needs to consider the effects of redshift-space distortions (RSD) and galaxy bias.First, the smooth component of the power spectrum takes the following form: where the factor B models the galaxy bias b g and the power spectrum amplitude variation.The term (1 +βµ 2 R(k)) 2 describes the effects of RSD at large scales [35].We have introduced the parameter β = f /b g , where f is the linear growth rate.The BAO oscillations are damped at high k in the power spectrum due to non-linear gravitational collapse.The non-linear evolution also causes subpercent shift in the BAO peak location [36].In order to revert the effects of non-linear evolution and thus enhance the statistical significance of the BAO signal, one can apply reconstruction, a technique to shift the density field back to its original positions [37].In that regard, the term R(k) in Eq. (3.9) that models the damping of non-linearities from RSDs on small scales can be modeled as [38] where Σ sm is a smoothing scale, and where the first line refers to the pre-reconstruction case and the second line to post-reconstruction 3 .The term F (k, µ, Σ FoG ) describes the Fingers of God effect [41] which is the elongation of observed structures in redshift space along the line of sight, primarily caused by the peculiar velocities on small scales.This term is defined as (see Eq. ( 27) in [42]) where Σ FoG is the streaming scale.It is a free parameter within the fitting template, and we assume Σ FoG = 10 Mpc/h in the fiducial cosmological model.The fitting template for the anisotropic matter power spectrum is finally where P shot ≡ 1/n and n is the mean number density of galaxies in the comoving volume.
The shot noise results from the discrete distribution of galaxies and it is calculated assuming a Poisson distribution of the galaxies.The Σ ⊥ and Σ ∥ parameters model the transverse and line-of-sight directions of the damping effects, respectively.Due to the lack of availability of direct measurements of µ in typical observations, the standard BAO analysis utilizes the Legendre multipoles of the anisotropic power spectrum which marginalize over µ.The monopole and the quadrupole moments of the fitting template of the anisotropic power spectrum take the following form: where L 2 is the Legendre polynomial of the second order.The expression ( accounts for the difference in the isotropic volume from the fiducial cosmology.Note that the factor ( r d,fid r d ) 3 is k-independent, and thus degenerate with (and can be subsumed in) the galaxy bias in the template-based analysis.The polynomial terms A l (k) are defined to be either one of these expressions where the first line refers to the pre-reconstruction case and the second line to post-reconstruction.The form of A l (k) for pre-and post-reconstruction is based on evaluating of the goodness of fit achieved by each term to the BOSS data [12].Since we are investigating models of modified gravity, it is not appropriate to use the power spectrum where the peak locations were enhanced using information from galaxy velocities (the "post-reconstruction power spectrum").This is because the reconstruction and its fitting template assume general relativity.Therefore, we study the monopole and quadrupole moments of the originally observed ("pre-reconstruction") power spectrum.In this case, there are 17 free parameters in the fitting template: a 0,1 , a 0,2 , a 0,3 , a 0,4 , a 0,5 , a 2,1 , a 2,2 , a 2,3 , a 2,4 , a 2,5 , In the fiducial cosmology, we choose the parameters of the template to take the following values: . Moreover, all the coefficients of A l (k) are set to 0 and, as mentioned before, the fiducial values of the alphas are unity by construction (α ∥ = α ⊥ = 1).We choose the fiducial value of Σ FoG following [22].The fiducial values for Σ ⊥ and Σ ∥ are motivated by [36].The value of b g depends on the survey and tracers under consideration and is presented in Table 2.The template and its fiducial parameters are used in the calculation of the covariance matrix, which will be discussed in the next section.

Fitting and extracting the BAO Signal
With the computed power spectrum in hand (Sec.3.2), and the description of how to model it (Sec.3.3), it is fairly straightforward to extract the BAO feature.We do so by fitting the BAO model that is described by parameters in Eq. (3.16).
First, we clarify that we adopt noiseless data.That is, we use the theoretically predicted power spectrum multipole moments, with error bars as described below but without adding stochastic noise.We adopt the noiseless data to reduce the sample variance in our results.To be clear, we do expect additional statistical error in the case with real data, but this statistical stochasticity operates independently of the biases caused by the insufficiently flexible template.It is precisely these latter effects that we wish to isolate and study.
To perform the fit, we need to define our likelihood.It takes the following form where P th is the concatenated vector of k-values of monopole and quadrupole moments of the power spectrum in the fitting template, and P data is that of the mock power spectra.Next, we need to specify C, the covariance matrix of multipoles.We start with the matter power spectrum of a small bin size of k and µ, which can be approximated as where ∆k is the width of k bins and ∆µ is width of the µ bin.The error bar here includes both cosmic variance and shot noise; note that the latter is implicitly included given that it appears in the expression for the anisotropic power spectrum in Eq. (3.12).The effective volume V eff is related to the measured physical volume by the equation: where n(z) is the number density of galaxies.See Sec. 4 for the chosen values of k and V eff for the test models.Assuming Gaussianity, the multipole covariance is [43] This implicitly includes the sub-covariance matrices for each multipole and those between different multipoles.We used a fixed covariance matrix given a galaxy survey, evaluated Table 2: Cosmological tracers in the adopted surveys.Here z has information about the redshift slice adopted, Ω Survey is the solid angle, V eff is the effective volume, n g is the average galaxy density, and b g is the galaxy bias that we adopted.

DESI SURVEY Tracer
z in the fiducial ΛCDM model; this is likely to be sufficiently accurate and also reflects the procedure adopted in typical compressed BAO analyses.We perform a global fit for all 17 parameters and effectively marginalize over 15 of them in order to obtain the posterior in the (α ∥ , α ⊥ ) plane.This approach stands in contrast to some other analyses that minimize over the other template parameters to constrain the alphas (e.g., [44]).While the two approaches appear to give comparable results in practice for ΛCDM model, the marginalization that we adopt is likely to be more robust when a wider range of cosmological models is considered.
To constrain the parameters p template , we employ the Markov chain Monte Carlo (MCMC) algorithm (emcee, [45]), ensuring convergence by adhering to the Gelman-Rubin convergence criteria.Specifically, we set a threshold for R values at less than 1.001 for each parameter.The MCMC algorithm uses the likelihood function as in Eq. (3.17).We adopt flat priors on each of the free parameters in the analysis.The α ∥ and α ⊥ parameters are varied in the range [0.8, 1.2], while the rest of the parameters were assigned flat priors considerably wider than their final posterior values.

Robustness Tests
Summarizing the procedure as outlined at the beginning of Sec. 3, we proceed as follows.For a given choice of Horndeski and background ΛCDM cosmological parameters, we first generate the power spectrum using EFTCAMB and compute its multipoles.Next, we fit these modified-gravity power spectra multipoles to the fitting template from Eqs. (3.13) and (3.14).The fitting template is computed based on a power spectrum assuming fiducial cosmology (section 3.4) and has seventeen free parameters that are listed in Eq. (3.16).Finally, we constrain all of these free parameters using the MCMC sampler.We repeat the procedure for 100 randomly sampled Horndeski models for each survey configuration described below.Furthermore, we investigate modified-gravity models where both background and perturbations significantly depart from ΛCDM.
We simulate data of the quality expected from the Stage-IV experiment Dark Energy Spectroscopic Instrument (DESI), its planned extensions DESI-II, and the Stage-V experiment modeled on the proposed MegaMapper survey.For DESI, we have initially considered three different traces that have been commonly used in recent BAO analyses: luminous red galaxy (LRG), emission line galaxy (ELG), and quasars (QSO).We found that the QSO constraints are relatively weak given their lower number density, so we only carried out the DESI simulation with the LRG and ELG.We estimated the properties of these tracers based on the in [46,47]. 4, and the tracer specifications are summarized in the top part of Table 2. Here, we quote the mean redshift, the redshift bin width (equal to 0.3 for both DESI tracers considered), as well as the solid angle, effective volume, number density, and galaxy bias that we assumed.In order to avoid the complexities of combining measurements from different redshift bins, for each tracer we only consider a single redshift bin.
In Table 2 we also show our specifications for DESI II and Megamapper.For these Table 3: Summary of the statistics of ∆α ∥ and ∆α ⊥ -differences between the measured and true values of the BAO peak-location parameters -based on 100 MCMC compressed analyses for different surveys and tracers.We show the mean and rms dispersion of each ∆α, along with the ratio of the rms dispersion in ∆α to the statistical error in the corresponding α, for each tracer type at specified redshifts in the DESI, DESI II, and MegaMapper surveys.We have carried out a battery of tests to validate our approach and code.First, we were able to reproduce the results of [22], where a very similar approach was adopted to study the robustness of several beyond-ΛCDM models to the compressed BAO analysis.We also studied the impact of different power-spectrum-smoothing methods to compute the power spectrum multiples.In particular, we compared the methodology utilized in the Barry code [51] 5 to that in [22].The two methods agree very well; see Appendix B for more details.We chose the configuration-space smoothing method in [22] as our fiducial approach.
Additionally, we tested modified-gravity models that simultaneously vary both background and perturbations.To achieve simultaneous variations in the background and perturbations, we allowed freedom in the functions Ω MG (t) and γ i (t).We specified these two functions according to power-law expressions in Eqs.(2.4)-(2.5),with corresponding amplitude and power-law coefficients sampled randomly from priors listed in Table 1.We further set the function Λ(t) to zero.The background evolution in our modified-gravity models is then specified by Ω MG (t) (see Eq. ( 4) in [30]), while the evolution of perturbations is controlled by it as well as by γ i (t).This framework allows us to explore the impacts on the compressed BAO analysis of modified-gravity models that significantly deviate from ΛCDM.The study of models where both background and perturbations vary from ΛCDM is not a primary goal of this paper, and these preliminary results are given in Appendix C.
Note there are two caveats for creating the mock power spectrum.First, the α ⊥ and α ∥ parameters in the mock power spectrum should be replaced by q ⊥ = D A (z) D A,fid (z) and q ∥ = H fid (z) H(z) .The reason is that the power spectrum has already been generated with the target true r d , so there is no need to further scale the r d from its fiducial value.Second, the units for the power spectrum P (k, µ) data need to include the same scaled Hubble constant h as that used in the fitting template (that is, the fiducial value of h adopted in the template).Only then can the correct definitions for α ⊥ and α ∥ be recovered from the fitting template.This is because there are no free parameters available to rescale h between the units of the mock power spectrum and those employed in the template.

Results
We now present our main results.We are interested in the constraints on α ∥ and α ⊥ , marginalized over the other 15 template parameters.Note that we are not particularly interested in the statistical error of the α parameters itself, given that it is dependent on our choice of the redshift bin along with all other survey specifications (which may end up being different in reality from what we assume here).Rather, we focus on the biases in the best-fit value of the alphas relative to the size of their statistical error.
In Fig. 3 we present the α ∥ −α ⊥ log-likelihood contour for one example randomly chosen modified-gravity model6 and for several survey configurations.The blue star indicates the true values of the alphas (which are quite close to unity for this particular Horndeski model), the red circle shows the best-fit value from our procedure described in Sec. 3, and the grey contours show the 68% and 95% credible contours for our fit, marginalized over all of the fit parameters.This figure already previews our key result, which is that the recovered alphas are in excellent agreement with their true values.In other words, the bias in the recovered alphas relative to their statistical values is well below 1 σ.This is true for both tracers of DESI, as well as DESI II and MegaMapper (the latter of which has a small forecasted error in the alphas even in a single redshift bin that we assumed).We therefore infer, from just one model for the moment, that the compression appears to be robust when applied to Horndeski models.
We next show the full statistics of the recovered α ∥ and α ⊥ parameters, applied to a sample of approximately 100 Horndeski models7 .We randomly sample these models from the Horndeski parameters ranges given in Table 1 and ΛCDM parameter values shown in Fig. 1.
To study the statistics of the BAO analysis robustness across the sample of models in our analysis, we define the shifts of the alpha parameters relative to their true values, We focus on the typical values of these shifts relative to typical statistical (measurement) errors in the corresponding alphas.
Table 3 presents the statistics of the shifts in the two alpha parameters derived from realizations conducted on various tracers and galaxy surveys.The mean and standard deviation, ∆α ∥ and σ(∆α ∥ ) (and same for α ⊥ ), are both found to be very small -between 10 −3 and 10 −2 .Perhaps more usefully, we also show the ratio of the typical deviations in the systematic biases in the alphas relative to their measurement errors, σ(∆α ∥ )/σ(α ∥ ) (and Figure 4: Distribution of the values of ∆α ∥ and ∆α ⊥ (red points) -differences between the measured and true values of the alphas-for mock Horndeski power spectra.The four panels refer to different tracers and survey configurations that we studied, while the black contour in each instance indicates the typical 68% and 95% credible contours in the α ∥ -α ⊥ plane.This plot illustrates that the biases ∆α ∥ and ∆α ⊥ are well within the typical statistical errors, reflecting the fact that the standard BAO analysis is robust against data assuming modified gravity theory.same for α ⊥ ).The typical values of the biases in the best-fit values of the recovered alphas are between 5% and 8% of their typical statistical errors.
We further illustrate our findings in Fig. 4, where the red points show the distribution of the biases in the alphas, ∆α ∥ and ∆α ⊥ , for the same four survey/tracer/redshift choices as in Fig. 3 and Table 2.In each instance, we also show the typical 68% and 95% credible contours in the corresponding alphas as the black contour.Given that the fixed fiducial P (k) is used to calculate the covariance matrix, the credible contours look similar across the same survey.Thus, typical contours suffice our illustrative purposes here.We again see that the recovered alphas are close to the true value.The recovered alphas remain well within the statistical error even in the cases with the most extreme biases.
Furthermore, we tested (in Appendix C) BAO fits to modified-gravity models where both background and perturbations depart from ΛCDM.In this preliminary analysis, we selected a handful of such models that range from being relatively close to agreeing with Planck data, to differing with it at >5σ.We found that the constraints are largely unbiased for these background-and-perturbations-varying models, and most true alphas lie within the fits' 95% credible contours.Furthermore, the biases in the recovered alpha parameters are small (<1σ) when base cosmological parameters are close (2-4σ away) to the fiducial cosmology.These results illustrate the robustness of the compressed BAO analysis even in modifiedgravity models with considerable differences between their background evolution and the fiducial cosmology of the template.However, we also find that when both background and perturbations are allowed to vary, the alpha parameters constraints are significantly weaker (∼ 1.5 − 2x larger contour sizes compared to models where only the perturbations are non-ΛCDM), and some models have notably poor constraints on the alpha parameters.While the additional biases when the background is varied may appear to be of concern, it is also the case that variations in the background are already limited and will be getting more so in the future, as the cosmological constraints on the geometrical measures that govern the background (from BAO, type Ia supernova, and other data) get better.We defer for the future a detailed study of compressed BAO fits in modified-gravity models where both background and perturbations are allowed to be significantly different from ΛCDM.
Overall, we have found that a compressed analysis of the BAO in Horndeski models returns accurate results for the key parameters α ∥ and α ⊥ , with biases that are well below the statistical error even for a future survey such as MegaMapper.

Conclusions and Discussion
Standard baryon acoustic oscillation (BAO) analyses compress the clustering information to the location of the BAO peak in the parameters α that are defined separately for the directions perpendicular and parallel to the line of sight and in each redshift bin.These compressed analyses utilize a physically motivated template to isolate the alphas from other information in the 3D power spectrum; the template is pre-computed and typically assumes a fixed cosmological model (e.g., ΛCDM with concordance values of cosmological parameters).It has been a long-standing question of just how robust this type of analysis is when considering more complex cosmological models.The robustness of this methodology has been tested for some specific departures from ΛCDM (see the Introduction), but not for a broad class of modifiedgravity models.There is some urgency to address this question since a principal goal of the ongoing Stage IV and forthcoming Stage V surveys is precisely to constrain modified gravity with BAO and RSD.
In this paper, we found that the compressed analysis is robust to a broad Horndeski class of modified-gravity models.Specifically, we have studied models where the perturbations are determined by Horndeski models, with eight additional free parameters that can freely vary, while the background is given by ΛCDM (with six standard cosmological parameters being varied).We have made use of the EFTCAMB implementation of Horndeski models to carry out the theoretical predictions, and have implemented our own analysis pipeline that follows the standard compressed-analysis approach.For each survey configuration, we studied 100 cosmological models that are not obviously ruled out (that is, that are in ≲5σ tension with Planck 2018 angular power spectrum data).Our results indicate that the biases in the recovered alphas are less than 10% of the statistical errors even for a Stage-V survey such as MegaMapper.Moreover, we have extended our study by considering (in Appendix C) modified-gravity models with both background and perturbations different significantly from ΛCDM.Even in that scenario, we found no significant biases (≲2σ) in the recovered alphas.However, this analysis was preliminary and will require further investigation to obtain precise statistical quantification of the fits for these models.
Overall, our findings, combined with previous work on other beyond-ΛCDM models (notably [22]), indicate that a compressed analysis based on a ΛCDM template remains remarkably robust with respect to the choice of the underlying cosmological model.
While we have established the robustness of the compressed analysis of a broad class of modified-gravity models, we did not cover all potential modifications of gravity (see [53] for a review).For example, one could further study beyond-Horndeski scalar-tensor models which have two additional free functions [54], or degenerate higher-order scalar-tensor (DHOST) theories which have higher-order equations of motion [55].One may also be interested in investigating modified-gravity models beyond scalar-tensor theory, such as higher-dimension, tensor-tensor, or tensor-vector-scalar theories [56].Current observational data have not shown statistically significant departures that would favor these models [57][58][59].However, forthcoming galaxy surveys such as DESI, LSST, Euclid and Roman telescopes, and the Stage-V spectroscopic instrument will provide significant improvement in statistical constraints that will make the observational analysis of these models more compelling.
Finally, we note that we have only tested the robustness of the compressed BAO analysis in this paper.This approach is fairly standard and well-established, but there now exist several more general methods that attempt to extract broadband information in the power spectrum beyond the BAO peak locations.These methods include BAO+RSD (i.e.fitting for f(z)σ 8 (z), e.g., [60]), ShapeFit [17], as well as direct modeling of the whole power spectrum [13][14][15][16].These methods are more general than the analysis that simply works off of α ∥ and α ⊥ and thus offer a greater potential to extract information from high-quality spectroscopic observations.These methods also make use of the broadband power spectrum which, as we have seen (e.g., in our Appendix A), is strongly impacted by modified gravity.There is therefore some level of urgency to study the robustness of these more ambitious methods to the underlying cosmological model.
In conclusion, while comparing the performance of different BAO and RSD methods on a wide range of cosmological models remains a priority, we have shown in this paper that the longest-established of such analyses -the compressed BAO analysis that uses a fixed template -is very robust in a wide range of modified-gravity cosmological models.
to the smoothed portion of each corresponding spectrum, O lin = P (k)/P smoooth (k).Note that varying γ 1,0 and s 1 parameters have minimal impact on the BAO feature compared to varying the background.This outcome aligns with expectations since the contributions from γ 1 (a) are relatively small when contrasted with observable uncertainties [32].Thus, we observe that P (k)/P smoooth, fid (k) is almost unchanged for these two parameters.In particular, we see that varying Ω MG,0 , s 0 , s 2 , and s 3 mostly affects the shape of the BAO feature.
We then compare the best-fitted parameters from the analysis applied to one randomly chosen Horndeski model with that applied to the fiducial ΛCDM model; the results are shown in Table 4.These two models have different input background cosmology (that is, the ΛCDM parameters assumed are different for the two models).The parameters are determined by adopting our analysis procedure on noiseless mock data, and using the least-squares fit implemented in the iminuit tool [61], assuming MegaMapper survey with LBG tracers at redshift of 2.25.The central column lists the best-fitted parameters for the analysis of the Horndeski model, while the right column details the parameters for the analysis of ΛCDM.Note the statistically significant deviation of the parameters α ∥ and α ⊥ in Horndeski analysis relative to the fiducial values of unity.This deviation is expected as the background cosmological parameters for the Horndeski model power spectrum are different from the fiducial-cosmology values, and lead to shifts in the BAO peak positions.
One other feature of note seen in Table 4 is a significant variation in the parameters that model the broadband power spectrum for the Horndeski analysis relative to those in the ΛCDM analysis.This is especially true for polynomial terms A l (k) (particularly higherorder terms), for Σ ⊥ and Σ ∥ which model the damping of BAO, and for the galaxy bias b g 2 .These differences between the two analyses are expected, and confirm that the amplitude characteristics of the power spectrum in Horndeski models are markedly distinct from those in ΛCDM.
Finally, we mention an additional caveat: modified-gravity models may predict different large-scale bulk flows which may in turn affect the BAO.These nonlinear effects manifest as a change in the amplitude of the BAO wiggles and the shape of the BAO feature with respect to those in the ΛCDM models but are not expected to affect the BAO positions.Such effects are modeled by means of IR resummation with time-sliced perturbation theory [62] in the direct modeling approach [13], and this modeling is expected to be accurate if higher-order calculations are included.However, the compressed analysis, which employs an exponential suppression term to model the BAO wiggles amplitude (as seen in Eq. 3.12), may not be sufficiently accurate to model the bulk-flow effects in modified gravity models.Our results indicate that adopting the exponential-suppression term in modified gravity models remains sufficiently accurate, as it does not bias the estimation of α ∥ and α ⊥ .Nonetheless, determining accurate analytical expressions for the exponential suppression term in power spectra under modified gravity models needs future investigation.

B Comparing Power Spectrum Smoothing Methods: Direct Interpolation vs. Indirect Approaches
Here we investigate the differences between the methods used to extract smooth components of the power spectrum in Barry code [51] and those adopted by [22].Barry employs direct interpolation on the power spectrum to achieve smoothing; specifically, it employs polynomial functions to dewiggle the power spectrum.On the other hand, the [22] approach involves a conversion to the configuration space for smoothing, followed by reconversion to the Fourier space; this procedure also smoothes the power spectrum.Fig. 7 visually contrasts these two methodologies.In the figure, we define the difference between the power spectra, ∆P (k), as: Here, P sm, config (k) denotes the smoothed power spectrum obtained from the configurationspace method presented in [22], while P sm, barry (k) corresponds to the spectrum derived using the Barry code.
To test the effects of these two smoothing methods, we apply them to the BAO compressed analysis for a single Horndeski model, with all other choices (e.g., survey specifications) being the same.The smoothing method in the Barry code reports α ∥ = 0.986 ± 0.048 and α ⊥ = 0.974 ± 0.010.In comparison, the smoothing method in the configuration space reports α ∥ = 0.985 ± 0.028 and α ⊥ = 0.975 ± 0.010.The differences between the best-fitted α ∥ and α ⊥ are less than 0.1% for these two smoothing methods, well below the statistical errors in the alphas.While the resulting peak locations agree extremely well between the two methods, we also see in Fig. 7 that the resulting amplitude of P (k)/P (k) smooth , which is of The matter power spectrum in the Horndeski models with varying γ 1,0 , γ 2,0 , γ 3,0 , s 1 , s 2 , and s 3 parameter, which are controlling perturbation evolution in these models.These plots illustrate the modifications in the shapes of the BAO signal by Horndeski models, similar to Fig. 5. less importance of the BAO analysis, also agrees well, with a typical difference around 0.5% between the two smoothing methods.Therefore, we have shown that the Barry and configuration-space smoothing methods agree very well, and the choice of which one to pick will have a negligible effect on the final results.
C Exploring modified-gravity models significantly deviating from ΛCDM in both background and perturbations Here we perform tests of the compressed BAO analysis on modified-gravity models in which both background and perturbations are allowed to depart from ΛCDM.We performed the compressed BAO analysis on ten models randomly selected according to the procedure in Section 3.5, so that the models are in the <5σ deviation from Planck's best-fit cosmology in the multi-dimensional parameter space.Two of the model fits had very poor constraints on the alpha parameters, and we omitted them from further analysis 10 .Moreover, we selected five more "extreme" models that are 5-6σ away from Planck's best fit; of these, three models had good fits on the alphas and we consider them further.For these preliminary tests of models where both the background and perturbations differ from ΛCDM, we only considered the MegaMapper survey settings.This plot shows the ratio of the power spectrum, P (k), to its smoothed versions using the smoothing method used in Barry framework [51] (solid line) and the indirect smoothing method (dashed line).
Fig. 8 shows the fits to the alphas for these 11 modified-gravity models.The eight models with <5σ deviation from Planck are shown as blue points, while the three models with 5-6σ deviation are shown as red points.We find that biases in the recovered alpha parameters are all below 2σ, and in 9 out of 11 cases, they are less than 1σ.The performance of the compressed BAO analysis therefore remains excellent even for these modified-gravity models with more freedom, at least based on this preliminary analysis with limited statistics.
To accommodate potential cases where the underlying model significantly deviates from ΛCDM, the performance of the BAO compressed analysis can also be improved by undertaking additional steps (not adopted in this paper).For example, one can extract information from BAO without relying on a precomputed template (e.g., [14,63]).In this case, at each step of the analysis, a linear power spectrum and high-order loop corrections are calculated to infer parameters of interest (instead of using a fixed template).Alternatively, one can add new degrees of freedom to the fiducial cosmology in the template to account for the possibility of modified gravity (by, for example, adding extra nuisance parameters to capture the impact of Ω MG (t) and γ i (a)).Investigating these possibilities will enable further overall robustness of the BAO methodology.
Figure 8: Distribution of the values of ∆α ∥ and ∆α ⊥ , representing differences between the best-fit values of the alphas, for a total of 11 randomly selected modified-gravity models where both background and perturbations differ from ΛCDM.The eight blue points are modifiedgravity models that are <5σ away from Planck's best-fit cosmology, while the three red points represent "extreme" models that deviate from Planck by 5-6σ.The black contours indicate the typical 68% and 95% credible contours in the α ∥ -α ⊥ plane.Note that 9 out of 11 models have biases in the alphas that are less than 1σ.We assume the MegaMapper survey settings here.

Figure 1 :
Figure 1: The parameters in the ΛCDM models for generating the mock power spectra data.The cosmological parameters were sampled within the 5σ credible interval of the Planck posterior.These parameters are sampled simultaneously with the EFTCAMB Horndeski model parameters.The red contour shows the 5σ credible (99.99994%) ellipses for the Planck data with TT+TE+EE+lowE constraint for comparison.

Figure 3 :
Figure 3: Likelihood contours in the α ⊥ − α ∥ plane with 68% and 95% confidence levels.The red and blue dots indicate the maximum-posterior and true values, respectively.The plot shows an analysis of a single cosmological model using the compressed analysis (for several survey choices).For our results reported in Sec. 4, we repeat this procedure about 100 times for different cosmological models.We consider four distinct galaxy surveys and tracers: the DESI survey with LRG tracers at z = 0.55 and ELG tracers at z = 1.25, the DESI II survey with LBG tracers at z = 2.25, and the MegaMapper Survey with LBG tracers at z = 2.25.

Figure 5 :
Figure5: Matter power spectrum in Horndeski models with varying Ω MG,0 (left column) and s 0 (right column).These two parameters control the background evolution of Horndeski models.In each column, the upper panel shows the power spectra divided by fiducial (fixed) power spectrum P (k)/P fid (k), the middle panel shows the power spectra divided by smoothed component of the fiducial power spectrum, P (k)/P smoooth, fid (k), while the bottom panel shows the spectrum divided by its own smooth components, O lin = P (k)/P smoooth (k).The lower panels therefore characterize the BAO feature; we see that the locations of the BAO do not change when the Horndeski model parameters are varied, but they do change once we also vary the base cosmological-model parameters.

3 Figure 6 :
Figure6: The matter power spectrum in the Horndeski models with varying γ 1,0 , γ 2,0 , γ 3,0 , s 1 , s 2 , and s 3 parameter, which are controlling perturbation evolution in these models.These plots illustrate the modifications in the shapes of the BAO signal by Horndeski models, similar to Fig.5.

Figure 7 :
Figure7: Comparison of power spectra with different smoothing methods.This plot shows the ratio of the power spectrum, P (k), to its smoothed versions using the smoothing method used in Barry framework[51] (solid line) and the indirect smoothing method (dashed line).

Table 1 :
The sampled ranges of Horndeski-model parameters that we considered to generate the mock power spectrum data.We assume a flat prior in all cases.

Table 4 :
Comparison of BAO compressed analysis best-fit parameters in the fiducial ΛCDM model and in a Horndeski model.Here, the set of parameters is for the anisotropic power spectrum fitting template.