Combined tracer analysis for DESI 2024 BAO

The Abacus simulations are high-resolution gravity-only N-body simulations [62]. We use a set of 25 simulation boxes from the AbacusSummit suite [63], each with a volume of $\rm(2h^{-1}Gpc)^{3}$ and $6912^{3}$ particles. These AbacusSummit simulations assumed the Planck 2018 $\Lambda$CDM cosmology, specifically the mean estimates of the Planck TT,TE,EE+lowE+lensing likelihood chains: $\Omega_{c}h^{2}=0.1200$, $\Omega_{b}h^{2}=0.02237$, $\sigma_{8}=0.811355$, $n_{s}=0.9649$, $h=0.6736$, and $w=-1$ [64]. This is also the fiducial cosmology used throughout the DESI DR1 analysis.

The Abacus-2 halos⁴⁴4Identified with the CompaSO halo finder [65]. are populated with galaxies using a flexible halo occupation distribution model (HOD) [66] which has been fitted to the galaxy two-point correlation functions of the final Early Data Release (EDR) [67, 68, 69] that included correction for all the systematics and included a detailed model for DESI focal plane effects. More details about the production and utilization of the mocks are provided in [39]. From these simulations, we create mock catalogs designed to mimic the survey realism, called “cutsky” Abacus-2 mocks, after applying the DESI footprint and survey selection function, including realistic redshift failures and targeting masks as described in [51]. As a caveat, the LRG DR1 mocks are produced using the simulation output at $z=0.8$, while the ELG DR1 mocks were produced using the simulation output at $z=0.950$ [51]. In comparison, the actual effective redshifts of the two DR1 samples (LRG3, ELG1) are very close: 0.92 and 0.95, respectively (Table˜1). Despite the offset in the output redshifts of the mocks, the cross-correlations between the two tracers reasonably agree between the mocks and the data on the BAO scale, as will be shown in Section˜5 (Fig.˜4). We, therefore, ignore the effect of this redshift offset.

As the final step, fiber assignment in the mocks is modeled using the DESI fiber assignment pipeline, which replicates the actual survey tiling, priority, and collision logic. It is important to simulate this effect as it might bias our observation and so our BAO measurements. For details, see [52, 70]. This is called the ‘altmtl’ (the Alternate Merged Target Ledgers [71]) Abacus DR1 mocks. Even though they are more computationally expensive, they are the most realistic simulations of the DR1 data as they accurately reproduce DESI fiber assignments. Fig.˜2 shows the clustering amplitude of the DR1 mocks for LRG, ELG, their cross correlation LRG$\times$ELG and their combination LRG$+$ELG.

The second type of mocks is based on the effective Zel’dovich approximation mocks [72], which we call EZmocks. EZmocks are almost as accurate as ${\tt Abacus}$ mocks on linear scales but less computationally expensive. We use 1000 realizations of DR1 EZmocks that are calibrated to produce two-point clustering statistics of Abacus-2 DR1, while constructed from a cubic volume of ${\rm(6h^{-1}Gpc)^{3}}$, as described in [39]. When applying the survey realisms to EZmocks, instead of applying and simulating fiberassign, we use the ‘fast-fiberassign (FFA)’ method that emulates the fiber assignment process by generating and averaging realizations from the average targeting probability as a function of number of overlapping tiles and local (angular) clustering, learned from the data. This ‘FFA’ method approximately reproduces the fiber-collisions pairwise incompleteness and is much faster than ‘altmtl’. We use these DR1 EZmocks only for the tests of the gain in BAO precision (Section˜5.1.2), which requires a greater number of independent simulations and needs to avoid a potential box replication effect ⁵⁵5To reach the full survey volume, the simulation box must be periodically replicated. This artificial tiling can introduce spurious large-scale correlations or suppress variance on scales approaching or exceeding the box size. [51]. The rest of the mock tests in this paper are performed using Abacus-2 DR1 mocks.

Our analysis requires an assumption of the fiducial cosmology for transforming the observational coordinates to the comoving coordinates and calculating the two-point correlation function, as well as in the BAO fitting. Whereas we use the Planck 2018 $\Lambda$CDM cosmology as the DESI fiducial cosmology (Section˜3.1), the robustness against the choice of fiducial cosmology is extensively demonstrated in [73].

Throughout this paper, we follow the definition of $\alpha_{\parallel},\alpha_{\perp}$ from the DESI DR1 BAO analysis [39].:

where $r_{d}$ is the sound horizon at the drag epoch, $D_{A}(z)$ is the angular diameter distance, $H(z)$ is the Hubble parameter, and $\alpha_{\parallel}$ and $\alpha_{\perp}$ are the dilations along and across the line of sight, respectively. Quantities with the superscript “fid” are measured in the fiducial cosmology that we introduced above. Because $D_{A}(z)$ also depends on $H(z)$, the two $\alpha$ parameters are correlated at some level. It is useful in some cases to express the 2-D dilation in terms of the isotropic $\alpha_{\rm iso}$ and anisotropic $\alpha_{\rm AP}$ distortion parameters defined as

where $\alpha_{\rm iso}$ measures the spherically averaged BAO scale, and $\alpha_{\rm AP}$⁶⁶6AP refers to the Alcock-Paczynski effect [74] i.e anisotropic distortion of the BAO feature. captures the anisotropy in the BAO feature.

The size and location of the BAO feature have been demonstrated to be approximately tracer-independent, especially after reconstruction [75, 76]. Therefore, the most naive approach to combining different BAO tracers is to concatenate the catalogs. However, given the different galaxy biases and sampling noise of the two tracers, an optimal approach would be to weight each tracer according to its signal-to-noise ratio when combining the catalogs. In this paper, we propose and test an approximate optimal weight assuming a scale-independent and isotropic signal-to-noise ratio. Once combined, we treat this as a single catalog, where the information it provides should be equivalent to or better than (especially if there is a gain during reconstruction) the combined information from all three clustering measurements available: the two auto-clustering measurements of the individual tracers and their cross-clustering. The cross-covariance among these three clustering statistics is naturally accounted for by treating them as a single catalog (at least within the approximation we use to derive the weight for the combination).

The ultimate information we want to derive from the combined tracer is the underlying matter density field (and the BAO feature in it) and the corresponding displacement field. LRG and ELG trace the underlying matter density field in a biased way. We construct an approximate maximum likelihood estimator of the matter density field, ignoring the anisotropies due to redshift-space distortions and scale-dependence of bias:

where $b_{1}$ and $b_{2}$ refer to the linear galaxy bias parameters for LRGs and ELGs, ${n}_{1}$ and $n_{2}$ are the observed number densities of LRG and ELG at a given location/pixel and $\bar{n}_{1}$ and $\bar{n}_{2}$ are the average number densities expected at given location/pixel, which are typically traced with randoms. The derivation for Eq.˜4.5 is presented in Appendix˜A.

Based on the estimator, we choose to weight each tracer with its estimated linear galaxy bias before merging the two catalogs. Then, the overdensity field of this combined catalog can be written as:

and the effective bias of this combined field is:

When comparing Eqs.˜4.6 and 4.7 with Eq.˜4.5, one can show

For simplicity, we ignore the position/redshift-dependence on $b_{\rm eff}$ when deriving the displacement field during reconstruction (Section˜4.4). But we account for the redshift-dependence of $b_{\rm eff}$ when updating the FKP weight for the combined tracer (Section˜4.2.2). The procedure for constructing the combined catalog can be summarized as follows:

1. 1.

   Estimate $\bar{n}_{\rm eff}(z)$ of the combined galaxy density field, assuming that each tracer is weighted with its estimated linear galaxy bias, and compute the new FKP weight based on it (Section˜4.2.2). This can be done without altering the clustering catalogs.

2. 2.

   In the clustering catalogs, weight LRGs with $b_{1}$ and ELGs with $b_{2}$, i.e., multiply the existing clustering weight by the corresponding bias for each galaxy. Weight the random particles in the same way.

3. 3.

   Renormalize the randoms so that the data-to-randoms weighted sum ratio is the same for both tracers. See Section˜4.2.3.

4. 4.

   Concatenate the two galaxy catalogs with the new weights (with biases, the new FKP weight). Do the same to the randoms (with the additional renormalization).

In the following subsections, we present practical details of each step for interested readers.

Following the pipeline defined in [39, 77], the correlation functions are computed using the Landy-Szalay estimator [78]

with $\mu\in[-1,1]$ being the cosine of the angle between the galaxy pair and the line of sight, $D,R$ refer respectively to points from the data and random catalogs, and $s$ is the pair separation (here we use bins with a width of $4h^{-1}\text{Mpc}$). To compute the two-point correlation we use the publicly available code pycorr¹⁰¹⁰10https://github.com/cosmodesi/pycorr [79] based on the Corrfunc¹¹¹¹11https://github.com/manodeep/Corrfunc engine [80, 81]. For our analyses, we convert the output to Legendre multipoles (specifically the monopole ($l=0$) and the quadrupole ($l=2$)). To extract the maximum out of the clustering measurement, we weigh each galaxy by a combination of weights, as described in Section˜4.2 (Eq.˜4.14 for the combined tracer, analogous to Eq.˜4.9 for regular single tracers).

Since the data reduction produces catalogs for each of the galactic caps, we choose to compute the correlation functions for NGC and SGC separately. Accordingly, we perform the catalog combination separately for NGC and SGC. Then the combination is performed by summing the pair counts computed in each region independently. This allows us to compare the consistency of the two regions and test for possible systematics. In some cases, an insufficient number of randoms can degrade the computation of the two-point statistics. To avoid this, we concatenate multiple random catalogs (typically 18) so that the number of randoms is more than 50 times the number of data galaxies.

We apply the density field reconstruction technique [82] on the catalogs considered in this paper, single and combined tracers, to partially recover the BAO feature that has been degraded due to structure growth and redshift-space distortions (RSD). We follow the DESI DR1 default reconstruction setup in [39], which was determined based on a set of extensive tests in [77, 83].

To summarize, we use pyrecon,¹²¹²12https://github.com/cosmodesi/pyrecon a Python package developed by the DESI collaboration, and adopt IterativeFFTReconstruction (hereafter, ‘IFFT’) that implements the iterative procedure described in [84], with the RecSym convention. ¹³¹³13RecSym is a choice to recover the large-scale anisotropy due to redshift-space distortions in the process of reconstruction. The overdensity field is smoothed by a Gaussian kernel of width $15h^{-1}\text{Mpc}$ for all of the tracers we consider in this paper. For reconstruction, we use the linear bias in Table˜2 and the growth rate in Table˜1. The value of linear bias, $b_{\rm eff}=1.6$, is computed from Eq.˜4.11 where the formula is integrated over the redshift range.

As a caveat, whereas the default reconstruction setup of the individual tracers is performed using the displacement field constructed over the full redshift range for LRG and ELG, i.e., $0.4<z<1.1$ for LRG and $0.8<z<1.6$ for ELG, the reconstruction of the combined catalog is performed using the displacement field estimated from $0.8<z<1.1$, i.e., the range where the two tracers overlap, as stated in [39]. The latter choice was made to input bias and growth rate that are better representative of the overlapping sample to the reconstruction pipeline, rather than using an effective bias and the growth rate averaged over a wide range of redshift (i.e. $0.4<z<1.1$ or $0.8<z<1.6$). The downside of this choice is a potentially greater boundary effect at $z=0.8$ when estimating the displacement field across the boundaries. The left panel of Fig.˜3 shows the mean correlation functions of five LRG$+$ELG mocks over $0.8<z<1.1$ using the displacement field constructed from $0.4<z<1.1$ (green-dotted, the convention used for LRGs) and from $0.8<z<1.1$ (blue, our default choice for LRG+ELG). The figure shows that the choice of the reconstruction redshift range between the two options has little impact on the reconstructed correlation function. Therefore, the comparisons of LRG and the combined tracer in this paper are minimally impacted by the different reconstruction redshift ranges.

In the right panel, we also test the effect of the Gaussian smoothing kernel for reconstruction. Given the effective number density of the combined tracer, greater than LRG and ELG, we test a more aggressive kernel $10h^{-1}\text{Mpc}$ (dashed blue line): the gain of using a smaller smoothing kernel does not appear obvious in the BAO feature of the monopoles, which mainly determines the constraint on $\alpha_{\rm iso}$. While there is a slight difference in the quadrupole, based on the little visible impact on the monopole, we adopt $15h^{-1}\text{Mpc}$ that is the default for DR1 LRG and ELG samples.

For the cross-correlation LRG$\times$ELG of LRG and ELG, we cross-correlate after each of LRG and ELG is separately reconstructed. One could also use the cross-correlation of the two tracers after reconstructing using the common, combined displacement field, especially given the potentially noisy displacement field from ELG alone. For a more conservative systematic/consistency check, we, however, decided to cross-correlate the individually reconstructed two fields, which can be subject to more errors.

In this work, we used semi-analytical covariance matrices for the 2-point correlation function multipoles produced with the RascalC code¹⁴¹⁴14https://github.com/oliverphilcox/RascalC [85, 86, 87, 88, 89]. It computes the covariance matrix terms based on an empirical 2-point function (but not 3- and connected 4-point functions) and applies a shot-noise rescaling to mimic non-Gaussian effects. We refer the readers to [89] for further details on the methodology.

The covariance matrix for the combined tracer was produced similarly to other single tracers, using the combined data and random catalogs¹⁵¹⁵15Scripts are available at https://github.com/misharash/RascalC-scripts/tree/DESI2024/DESI/Y1/comb.. Therefore, the LRG, ELG and LRG$+$ELG RascalC covariance matrices in this paper are the same as in the main DESI 2024 galaxy BAO analysis [39]. The covariance matrix of the cross-correlation function (LRG$\times$ELG) was produced based on data correlation functions and random catalogs according to Appendix A of [89]¹⁶¹⁶16Scripts are available at https://github.com/misharash/RascalC-scripts/tree/DESI2024/DESI/Y1/cross., using the shot-noise rescaling values calibrated with jackknives on LRG and ELG auto-correlations. Table˜3 summarizes the covariances we use in this paper.

The BAO fitting procedure for the DESI collaboration is based on the extensive tests that are summarized in [38]. The model consists of a combination of a physically motivated model from quasi-linear theory and a parameterized model to marginalize over non-linearities that may otherwise affect our measurements of the BAO scale. For interested readers, [39, 38] explain the fitting model for different choices of the reconstruction convention, how the power spectrum multipoles are transformed into configuration space multipoles, the broadband modeling and its marginalization, etc. Here we simply show in Table˜4, a few of the parameters of interest.

To summarize, we use a mixture of flat ($\alpha_{\rm iso},\alpha_{\rm AP}$) and Gaussian priors (damping terms). To isolate the BAO feature, the range of galaxy separation is restricted to $48<s<152h^{-1}\text{Mpc}$. The fit is performed with the publicly available code on two-point correlation functions of combined galactic caps (NGC+SGC) using the desilike¹⁷¹⁷17https://github.com/cosmodesi/desilike. Fitting methods available are MCMC sampling (e.g. emcee [90], pocoMC [91]) or posterior profiling (MINUIT [92]). The software offers the possibility to fit in several parameter spaces. In this paper, we chose to present our results in $\alpha_{\rm iso}-\alpha_{\rm AP}$ or $\alpha_{\parallel}-\alpha_{\perp}$.

For the priors on the damping parameters, the central values denoted by the superscript ‘fid’ are tracer-specific and they are shown in Table˜2. These fiducial values are motivated by a combination of theoretical calculations, measurements of the cross-correlation between the initial and post-reconstruction density fields, and fits to mock catalogs. Finally, the total number of free parameters for our 1D and 2D fits is 7 and 13, respectively, for $\xi(r)$. The corresponding numbers of degrees of freedom for our 1D and 2D fits are 19 and 39.

We first apply our pipeline on the 25 Abacus-2 DR1 mocks and estimate the expected level of systematics and the gain from the combination of catalogs.

In Fig.˜4 we show the mean of the correlation function for the 25 realizations of pre- and post-reconstruction for LRGs, ELGs, LRG$+$ELG, and LRG$\times$ELG. The amplitude of clustering is proportional to the product of the galaxy biases of the two samples, and we can estimate the galaxy bias of LRG$+$ELG and LRG$\times$ELG from LRG ($b_{\rm LRG}=2$) and ELG ($b_{\rm ELG}=1.2$) ¹⁸¹⁸18The default bias values for LRG and ELG are taken from [39].. One can see a clear BAO feature around the 100 $h^{-1}\text{Mpc}$ with a slightly sharper peak after reconstruction (a sign that some of the non-linearities have been successfully corrected). This is especially the case for LRG and LRG$+$ELG. The ELG reconstruction was likely hindered by its low completeness, which should be improved with future data releases of the DESI survey.

The last step of our analysis is to apply the combined tracer analysis on the DR1 data and investigate the gain and potential systematics. In Table˜7, we present the results of BAO fits applied to the DESI DR1 data. First, Fig.˜8 repeats one of the sanity checks presented in [39]: pre- and post-reconstruction BAO measurements for LRG$+$ELG, comparing DESI DR1 measurements (stars) with the 25 Abacus-2 DR1. Therefore, the reconstruction efficiency of DESI DR1 LRG$+$ELG is consistent with what is expected from the mocks.

In terms of a gain in precision between LRG$+$ELG and LRGs, Table˜7 shows that, for $\alpha_{\rm iso}$ we find a gain of $5\%$ pre recon and $11\%$ post reconstruction, while the gain for $\alpha_{\rm AP}$ is $\sim 2.2\%$ pre recon and $\sim 7.0\%$ post reconstruction. Fig.˜9 visualizes that post-reconstruction BAO measurements between LRG and LRG$+$ELG (y-axis) as well as LRG and LRG$\times$ELG (x-axis) fall within the distribution of the mocks. Especially, the gain in precision of the combined tracer observed from the data (right panels) is consistent with the gain in terms of $\langle{\sigma^{L}_{\alpha}}\rangle$ of the Abacus-2 mocks.

Fig.˜7 shows differences in $\alpha$’s from pairs of the BAO tracers, including LRG$\times$ELG, comparing those of the mocks to the data. All differences from DR1 data are within the range of the histograms, therefore consistent with zero offsets, and therefore no indication of a tracer-dependent systematics/physics²⁴²⁴24Since we are testing a single realization, DR1, we compare the DR1 constraint with the entire width of the distribution, instead of the width divided by $\sqrt{(}25)$.. As a caveat, while LRG$\times$ELG returns consistent precisions compared to the mocks, the data returns $\chi^{2}$ much higher than the average $\chi^{2}$ of the mocks in Table˜5. Although this implies $\xi_{\rm LRGxELG}$ may need a more careful inspection in terms of the goodness of the fits, we defer such a test to future study.