The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: cosmological implications from multitracer BAO analysis with galaxies and voids

ABSTRACT

We construct cosmic void catalogues with the dive void finder upon SDSS BOSS DR12 and eBOSS DR16 galaxy samples with BAO reconstruction, and perform a joint BAO analysis using different types of galaxies and the corresponding voids. The BAO peak is evident for the galaxy–galaxy, galaxy–void, and void–void correlation functions of all data sets, including the ones cross-correlating LRG and ELG samples. Two multitracer BAO fitting schemes are tested, one combining the galaxy and void correlation functions with a weight applied to voids, and the other using a single BAO dilation parameter for all clustering measurements. Both methods produce consistent results with mock catalogues, and on average ∼10 per cent improvements of the BAO statistical uncertainties are observed for all samples, compared to the results from galaxies alone. By combining the clustering of galaxies and voids, the uncertainties of BAO measurements from the SDSS data are reduced by 5–15 per cent, yielding 0.9 per cent, 0.8 per cent, 1.1 per cent, 2.3 per cent, and 2.9 per cent constraints on the distance |$D_{_{\rm V}}(z)$|⁠, at effective redshifts 0.38, 0.51, 0.70, 0.77, and 0.85, respectively. When combined with BAO measurements from SDSS MGS, QSO, and Lyα samples, as well as the BBN results, we obtain |$H_0 = 67.58 \pm 0.91\, {\rm km}\, {\rm s}^{-1}\, {\rm Mpc}^{-1}$|⁠, Ω_m = 0.290 ± 0.015, and |$\Omega _\Lambda h^2 = 0.3241 \pm 0.0079$| in the flat-ΛCDM framework, where the 1σ uncertainties are around 6 per cent, 6 per cent, and 17 per cent smaller respectively, compared to constraints from the corresponding anisotropic BAO measurements without voids and LRG–ELG cross-correlations.

methods: data analysis, cosmological parameters, distance scale, large-scale structure of Universe

1 INTRODUCTION

The baryon acoustic oscillation (BAO) is a pattern in the clustering of matter, imprinted by sound waves in the primordial baryon–photo plasma of the Universe. It characterizes the sound horizon, a particular scale that can be used as a ‘standard ruler’ for measuring cosmological distances (Blake & Glazebrook 2003). In addition, it can be measured from both the cosmic microwave background (CMB) anisotropies (e.g. Bennett et al. 2003), as well as the distributions of large-scale structures at lower redshifts (Seo & Eisenstein 2003). The evolution of this scale across cosmic time traces the expansion history of the Universe, it is thus useful for constraining the properties of different cosmological components, such as the nature of dark energy (e.g. Weinberg et al. 2013; Aubourg et al. 2015).

Since the first clear detections of BAO from galaxy clustering achieved by the Two Degree Field Galaxy Redshift Survey (Cole et al. 2005) and the Sloan Digital Sky Survey (SDSS; Eisenstein et al. 2005), measurement of the BAO scale has been a key science goal for many large-scale spectroscopic surveys, such as the 6dF Galaxy Survey (6dFGS; Jones et al. 2004), the WiggleZ Dark Energy Survey (Drinkwater et al. 2010), the SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013), and SDSS-IV Extended BOSS (eBOSS; Dawson et al. 2016), the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016), as well as future surveys like the 4-metre Multi-Object Spectroscopic Telescope (4MOST; de Jong et al. 2019) and Euclid (Laureijs et al. 2011).

Nowadays, the precision of BAO measurements has achieved one per cent level with the distributions of luminous red galaxies (LRGs; Beutler et al. 2017; Ross et al. 2017), Moreover, the BAO feature has been detected in the clustering of different types of tracers, including star-forming emission-line galaxies (ELGs; Raichoor et al. 2021), quasi-stellar objects (QSOs; Ata et al. 2018), Lyman α (Lyα) forests (Busca et al. 2013), and cosmic voids (Kitaura et al. 2016a). With the fast expansion of the spectroscopic databases, it is even possible to perform BAO measurements in a multitracer manner, with the cross-correlations between different tracers taken into account, including but not limited to LRGs and ELGs (Wang et al. 2020; Zhao et al. 2021b), Lyα forests, and QSOs (du Mas des Bourboux et al. 2020), as well as galaxies and voids (Nadathur et al. 2019b; Zhao et al. 2020). In the multitracer approach, it has been shown that statistical uncertainties of various cosmological constraints due to cosmic variance effects are significantly reduced (see also McDonald & Seljak 2009; Seljak 2009; Blake et al. 2013; Favole, Sapone & Silva Lafaurie 2019).

Multitracer analysis typically requires observations of different types of tracers in the same cosmic volume, which are not always possible due to the inhomogeneity of tracer abundances at different redshifts and limitations of instruments. However, as an indirect type of tracers, cosmic voids can in principle be identified from all existing matter tracer catalogues without carrying out additional observations. Thus, it is almost always possible to perform multitracer cosmological measurements with voids. Voids are large regions that are devoid of luminous objects, and trace underdensities of the matter density field (see van de Weygaert & Platen 2011, for a review), though there are various different void definitions and identification algorithms (e.g. Colberg et al. 2008; Zhao et al. 2016; Cautun et al. 2018, and references therein). Therefore, the distribution of voids encodes the information from troughs of the density field, which are complementary to the peaks traced by matter tracers, such as galaxies and QSOs. This is true even for a linearized matter tracer distribution, which can be obtained by reconstruction algorithms (e.g. Padmanabhan et al. 2012), as these objects always trace only density peaks. As the result, the combination of voids and matter tracers is foreseen to provide a more detailed description of density field and yield better cosmological constraints.

This has been confirmed by Zhao et al. (2020) in terms of BAO measurements, which makes use of the void definition based on Delaunay triangulation (Delaunay 1934) for capturing the BAO feature from underdensities, using the Delaunay trIangulation Void findEr¹ (dive; Zhao et al. 2016). These voids are allowed to overlap with each other, resulting in a high void abundance, which turns out to be crucial for evidencing the BAO peak in void clustering (Kitaura et al. 2016a). Actually, it is so far the only void definition that can be used for BAO constraints. By combining the two-point correlation functions (2PCFs) of galaxies and voids, Zhao et al. (2020) presents a 10 per cent improvement on the BAO measurement precision based on studies using simulations, which is also achieved with the BOSS LRG data at low redshifts.

In this study, we extend the work in Zhao et al. (2020), and aim at a joint BAO analysis with all BOSS and eBOSS galaxies, as well as the corresponding voids. The measurements are validated using previous studies with the same data sets, including those for BOSS LRGs (Alam et al. 2017), eBOSS LRGs (Gil-Marín et al. 2020; Bautista et al. 2021), and eBOSS ELGs (Tamone et al. 2020; de Mattia et al. 2021; Raichoor et al. 2021). In addition, measurements from the SDSS main galaxy sample (MGS; Ross et al. 2015), eBOSS QSOs (Neveux et al. 2020; Hou et al. 2021), and Lyα forests (du Mas des Bourboux et al. 2020) are also included for cosmological parameter constraints. A crucial component of these studies is the construction of approximate mock catalogues (Farr et al. 2020; Lin et al. 2020; Zhao et al. 2021a) and galaxy mocks based on N-body simulations for assessing theoretical systematic errors (Avila et al. 2020; Smith et al. 2020; Alam et al. 2021b; Rossi et al. 2021). Our results are then compared to the cosmological analysis with the same data sets, but not including voids (Alam et al. 2021a, hereafter eBOSS21). In practice, existing SDSS BAO measurements and the cosmological interpretations are taken from the SDSS website.²

The data sets used in this work are described in Section 2, followed by detailed explanations of the analysis methodology in Section 3, including the construction of void catalogues and clustering measurement algorithms, as well as BAO measurement and cosmological parameter constraint schemes. The method is then validated using mock catalogues in Section 4, and analysis with the SDSS data are presented in Section 5. Finally, the results are concluded in Section 6.

2 DATA

We introduce in this section the data we use for BAO measurements, i.e. the LRG and star-forming ELG samples from BOSS and eBOSS observations – which are the same samples as used in eBOSS21 – as well as the corresponding mock catalogues for covariance matrix estimations and methodology validations, including both approximate mocks and catalogues from N-body simulations.

2.1 SDSS-III BOSS DR12 LRGs

The BOSS (Dawson et al. 2013) of SDSS-III (Eisenstein et al. 2011) measured the spectra of over 1.3 million LRGs between 2009 and 2014, using double-armed spectrographs (Smee et al. 2013) as well as the 2.5-metre Sloan Telescope (Gunn et al. 2006) at the Apache Point Observatory. The complete data set of BOSS LRGs is released as part of the SDSS Data Release 12³ (DR12; Alam et al. 2015). This LRG sample consists of two major populations, LOWZ and CMASS, that are targeted using different algorithms (Reid et al. 2016), for galaxies with redshifts ≲ 0.4, and from 0.4 to 0.7, respectively. Moreover, LRGs with redshifts measured by SDSS-I/II (Abazajian et al. 2009) are also included in DR12, as the BOSS target selections are designed to extend that of previous SDSS LRG samples (Eisenstein et al. 2001).

The final BOSS DR12 LRG catalogue for large-scale structure analysis is constructed by combining subsamples with different target selection rules, as the differences in clustering amplitudes are small (Reid et al. 2016; Alam et al. 2017). The footprints of this LRG sample for both northern and southern galactic caps (NGC and SGC, respectively) are shown in Fig. 1, with a sky coverage of nearly 10 000 deg² in total. Observational systematics, including photometric and spectroscopic effects, are corrected by various weights, and the total systematic weight is given by (Reid et al. 2016; Ross et al. 2017)

$$\begin{eqnarray*} w_{\rm tot} = w_{\rm sys} (w_{\rm cp} + w_{\rm noz} - 1), \end{eqnarray*}$$

(1)

where w_sys indicates the total angular photometric weights, while w_cp and w_noz are for correcting the fibre collision effect and redshift failures, respectively. The overall weight for clustering measurements is then

$$\begin{eqnarray*} w_{\rm all} = w_{\rm tot} w_{_{\rm FKP}} , \end{eqnarray*}$$

(2)

with |$w_{_{\rm FKP}}$| being the FKP weight for reducing clustering variances (Feldman, Kaiser & Peacock 1994), i.e.

$$\begin{eqnarray*} w_{_{\rm FKP}} = \frac{1}{1 + n(z) P_0} . \end{eqnarray*}$$

(3)

Here, n(z) is the radial comoving number densities of galaxies (weighted by w_tot, also known as the radial selection function, see Fig. 2), and a P₀ value of |$10\,000\, h^{-3}\, {\rm Mpc}^3$| is used for BOSS LRGs (Reid et al. 2016).

Figure 1.

Sky coverage of BOSS DR12 and eBOSS DR16 galaxies, for both NGCs (upper) and SGCs (bottom).

Figure 2.

Weighted comoving number densities of BOSS DR12 and eBOSS DR16 galaxies, with all the photometric and spectroscopic systematic weights included. Shadowed areas indicate redshift ranges of galaxy samples used in this work. The comoving distances and volumes are computed in the flat ΛCDM cosmology with Ω_m = 0.31.

Following Alam et al. (2017) and eBOSS21, we divide the LRG catalogue into two redshift bins. The redshift ranges are 0.2 < z < 0.5 and 0.4 < z < 0.6, and the samples are labelled ‘LRG (a)’ and ‘LRG (b)’, respectively. The effective redshifts of the samples are defined as the weighted mean redshift of all galaxy pairs with separations in the range (s_low, s_high), as is done in eBOSS analysis (e.g. Bautista et al. 2021):

$$\begin{eqnarray*} z_{\rm eff} = \left. \sum _{{\substack{i,j \\ s_{\rm low} \lt s_{\mathrm{ ij}} \lt s_{\rm high}}}}^{N_{\rm gal}} w_{{\rm all}, i} \, w_{{\rm all}, j}\, \frac{z_i + z_j}{2} \middle / \quad \sum _{ {\substack{i,j \\ s_{\rm low} \lt s_{\mathrm{ ij}} \lt s_{\rm high}}}}^{N_{\rm gal}} w_{{\rm all}, i} \, w_{{\rm all}, j} \right. \, , \end{eqnarray*}$$

(4)

where N_gal denotes the total number of galaxies in the given redshift ranges, and s_ij indicates the separation between every galaxy pair. We choose a separation range of |$[50, 160]\, h^{-1}\, {\rm Mpc}{}$| for effective redshift evaluations of all our samples, as it is close to the range used for our BAO measurements (see Section 4.2). The numbers of tracers and effective redshifts of our LRG samples are listed in Table 1. Note that when keeping two significant figures, the effective redshifts turn out to be identical to the ones used in BOSS DR12 analysis (e.g. Alam et al. 2017), though their definitions are different. Furthermore, BOSS LRGs with redshifts above 0.6 are not treated as an individual sample, but combined with eBOSS LRGs, as is detailed in the next section.

Table 1.

Redshift ranges and effective redshifts of the SDSS galaxy samples used in this paper, as well as the numbers of galaxies for each galactic cap.

Sample	Redshift		Galactic
	Range	zeff	Cap	Ngal
LRG (a)(BOSS DR12)	(0.2, 0.5)	0.38	NGC	429182
			SGC	174819
			Total	604001
LRG (b) (BOSS DR12)	(0.4, 0.6)	0.51	NGC	500872
			SGC	185498
			Total	686370
LRG (c) (CMASS+eBOSS DR16)	(0.6, 1.0)	0.70	NGC	255741
			SGC	121717
			Total	377458
ELG (eBOSS DR16)	(0.6, 1.1)	0.85	NGC	83769
			SGC	89967
			Total	173736

2.2 SDSS-IV eBOSS DR16 LRGs and ELGs

As the extension of BOSS, the eBOSS (Dawson et al. 2016) of SDSS-IV (Blanton et al. 2017) probes galaxy distributions at higher redshifts. It measured the spectra of over 0.2 million LRGs and slightly more star-forming ELGs from 2014 to 2019, using the same instruments as BOSS. The eBOSS LRG and ELG targets are selected to overlap with each other, as well as with the BOSS CMASS LRGs, in both angular and radial directions (Prakash et al. 2016; Raichoor et al. 2017), thus permitting cross-correlation measurements. Moreover, the target selection algorithm for LRGs is designed to avoid BOSS CMASS targets (Prakash et al. 2016).

The final eBOSS large-scale structure catalogues are released as DR16⁴ (Ahumada et al. 2020). The redshift ranges of LRGs and ELGs are 0.6 < z < 1.0 and 0.6 < z < 1.1, with the sky coverages of around 4200 and 720 deg², respectively. The footprints are shown in Fig. 1. Similar to the case of BOSS data, observational systematics of eBOSS galaxies are corrected by weights accounting for angular photometric effects, fibre collisions, and redshift failures (Ross et al. 2020; Raichoor et al. 2021). However, the combination of spectroscopic weights are different from that of equation (1), and the total weight for eBOSS galaxies is defined as

$$\begin{eqnarray*} w_{\rm tot} = w_{\rm sys} w_{\rm cp} w_{\rm noz}. \end{eqnarray*}$$

(5)

As stated before, eBOSS LRGs are combined with BOSS CMASS galaxies with z > 0.6 for cosmological analysis (Ross et al. 2020). We denote this combined LRG sample by ‘LRG (c)’ hereafter. The weighted radial distributions of eBOSS DR16 galaxies are shown in Fig. 2. With P₀ for FKP weights (see equation (3)) being 10 000 and |$4000\, h^{-3}\, {\rm Mpc}^3$| for the ‘LRG (c)’ and eBOSS ELG samples (labelled ‘ELG’), respectively (Ross et al. 2020; Raichoor et al. 2021), the effective redshifts given by equations (2–5) are listed in Table 1. In addition, the effective redshift computed in the same way, for the cross-correlation between these two samples, is |$z_{\rm eff}^\times = 0.77$|⁠, which is consistent with the value reported in Wang et al. (2020).

The eBOSS DR16 data includes also 0.34 million QSOs (Ross et al. 2020), which can in principle be used for the multitracer BAO measurements with voids as well. However, since the comoving number density of the QSOs is relatively low (⁠|$\sim 2 \times 10^{-5}\, h^3\, {\rm Mpc}^{-3}$|⁠), the radii of voids identified from the QSO sample are peaked at around |$35\, h^{-1}\, {\rm Mpc}$|⁠. The exclusion effect of such large voids may contaminate the BAO signature (see Liang et al. 2016, and discussions in Section 3.2.3), which requires careful inspections to avoid biases of BAO measurements. We leave a thorough multitracer BAO study with QSO voids to a future work.

2.3 DR12 MultiDark-Patchy mock catalogues

We rely on the DR12 MultiDark-Patchy (MD-Patchy) mock catalogues⁵ (Kitaura et al. 2016b) for estimating the covariance matrices of clustering measurements from the BOSS DR12 data, and validating our BAO fitting method (see Section 3.4). These mocks are created based on halo catalogues generated at 10 different redshifts in the range of 0.2 < z < 0.75, that account for the evolution of clustering. The halo catalogues are constructed with the cosmological parameters listed in Table 2, using the PerturbAtion Theory Catalogue generator of Halo and galaxY distributions (patchy; Kitaura, Yepes & Prada 2014), followed by halo mass assignment using the Halo mAss Distribution ReconstructiON method (hadron; Zhao et al. 2015). Both steps require calibrations with the BigMultiDark simulation (BigMD; Klypin et al. 2016, see also Section 2.5). Then, galaxy catalogues are constructed using the SUrvey GenerAtoR code (sugar; Rodríguez-Torres et al. 2016), which applies survey geometry to the halo samples, with observational effects such as stellar mass incompleteness and fibre collisions taken into account.

Table 2.

Flat-ΛCDM cosmology models used in this paper. The ‘fiducial‘ cosmology is for coordinate conversions and BAO fits of the SDSS data, MD-Patchy, and EZmock catalogues. The rest of the cosmology models are the ones for the construction of the corresponding data sets. Moreover, the BAO fits for the BigMD and OuterRim simulations are performed with their own cosmologies. Here, r_d denotes the comoving sound horizon at the drag epoch, computed using camb (see also Section 3.4.1).

Parameter	Fiducial	BigMD and	EZmock	OuterRim
		MD-Patchy
h	0.676	0.6777	0.6777	0.71
Ωm	0.31	0.307 115	0.307 115	0.264 79
Ωbh2	0.022	0.022 14	0.022 14	0.022 58
σ8	0.8	0.8288	0.8225	0.8
ns	0.97	0.9611	0.9611	0.963
Σmν (eV)	0.06	0	0	0
|$r_{\rm d}\, ({\rm Mpc})$|	147.78	147.66	147.66	149.35

The patchy code implements accurate models of structure formation based on the Augmented Lagrangian Perturbation Theory (Kitaura & Heß 2013), as well as galaxy biases that encodes nonlinear, stochastic, and nonlocal effects (Kitaura et al. 2014). As the result, MD-Patchy mocks reproduce the clustering of the BOSS DR12 data down to the scale of a few |$h^{-1}\, {\rm Mpc}$|⁠, including three-point statistics (Kitaura et al. 2016b). There are 2048 independent realizations of DR12 MD-Patchy mocks in total, for each galactic cap. We use only 1000 of them (the indices are from 0001 to 1000) in this work, for both ‘LRG (a)’ and ‘LRG (b)’ samples.

2.4 DR16 EZmock catalogues

To assess the clustering covariances of the eBOSS DR16 data, including those of the cross-correlations between LRGs and ELGs, we use the multi-tracer EZmock catalogues⁶ (Zhao et al. 2021a) generated using the effective Zel’dovich approximation mock generator (EZmock; Chuang et al. 2015a). These mocks build upon the Zel’dovich approximation (Zel’dovich 1970) and effective descriptions of tracer biases that need to be calibrated with simulations or observations, which is similar to that of patchy. EZmock is highly efficient, and still achieves comparable precision as patchy on the same scales, for both two- and three-point clustering statistics (Chuang et al. 2015b).

The DR16 EZmock catalogues are calibrated directly with the auto correlations of CMASS (for the combined LRG sample) and eBOSS data in different redshift bins. In addition, the same white noises of the density fields are used for different tracers, to account for their cross-correlations and covariances, as they reside in the same cosmic volume. Mocks generated at different redshifts – for reproducing the evolution of clustering – are then combined and applied survey geometry using the make_survey⁷ (White, Tinker & McBride 2014) and brickmask⁸ (Zhao et al. 2021a) codes.

Furthermore, various photometric and spectroscopic systematic effects are applied to the mocks, with corrections given by weights defined in the same way as those of the eBOSS data (see equation (5)). The clustering measurements of the final EZmock catalogues are in good agreement with those of the corresponding data on scales down to a few |$h^{-1}\, {\rm Mpc}$|⁠, including cross-correlations in configuration space, and the discrepancies are generally less than 1σ (Zhao et al. 2021a). We use all the 1000 realizations of DR16 EZmock catalogues in this paper for the ‘LRG (c)’ and ‘ELG’ samples, as well as their cross-correlations.

2.5 N-body simulations

Since the nonlinear damping of the BAO peak (e.g. Crocce & Scoccimarro 2008; Prada et al. 2016) may not have been fully captured by the approximate mocks, which is a known issue for the DR16 EZmock catalogues (Bautista et al. 2021; Raichoor et al. 2021), we use accurate mocks from high-resolution N-body simulations to model the BAO damping and validate our fitting method (see Section 3.4.1). To this end, we rely on galaxy catalogues constructed using N-body simulations at similar redshifts as our data samples, which reproduce the clustering statistics of the data.

For the LRG samples, we use galaxy catalogues from the BigMD simulation. They are actually the reference catalogues for the calibrations of DR12 MD-Patchy mocks (see Section 2.3), and are fitted to the BOSS DR12 data using the subhalo abundance matching (SHAM; e.g. Nuza et al. 2013) algorithm (Rodríguez-Torres et al. 2016). The BigMD simulation contains 3840³ particles with the mass of |$2.36 \times 10^{10}\, h^{-1}\, {\rm M}_\odot$|⁠, in a box with the side length of |$2.5\, h^{-1}\, {\rm Gpc}$| (Klypin et al. 2016). It is run with the TreePM N-body code gadget-2 (Springel 2005), using the cosmological parameters listed in Table 2. The halo catalogues used for SHAM are generated by the rockstar algorithm (Behroozi, Wechsler & Wu 2013). We choose three of the 10 BigMD galaxy catalogues for MD-Patchy calibrations, as is shown in Table 3, whose redshifts are the closest to the effective redshifts of our SDSS LRG samples (see Table 1).

Table 3.

N-body simulation galaxy catalogue used in this work, as well as the SDSS data that they are compared to (see Table 1). Properties of the simulation catalogues are also listed, including the side lengths of the periodic boxes, and redshifts and comoving number densities of the catalogues.

Simulation	SDSS	Boxsize		Number density
Sample	Sample	(⁠|$h^{-1}\, {\rm Gpc}$|⁠)	Redshift	(⁠|$h^3\, {\rm Mpc}^{-3}$|⁠)
BigMD (a)	LRG (a)	2.5	0.392	2.66 × 10−4
BigMD (b)	LRG (b)	2.5	0.505	4.28 × 10−4
BigMD (c)	LRG (c)	2.5	0.638	1.98 × 10−4
OuterRim	ELG	3	0.865	1.74 × 10−3

We use a galaxy catalogue generated using the halo occupation distribution (HOD; e.g. Berlind et al. 2003; Zheng et al. 2005) model (‘HOD-1’ in Avila et al. 2020), based on the OuterRim simulation (Heitmann et al. 2019), for our ELG sample. The OuterRim simulation is performed using the Hardware/Hybrid Accelerated Cosmology Code (Habib et al. 2016), with 10 240³ particles in the volume of |$(3\, h^{-1}\, {\rm Gpc})^3$|⁠. The particle mass is |$1.89 \times 10^9\, h^{-1}\, {\rm M}_\odot$|⁠, and the cosmological parameters are shown in Table 2. Halo catalogues are constructed using the Friends-of-Friends (Davis et al. 1985) algorithm, and the one used in this work is at z = 0.865. The HOD catalogue is fitted to the eBOSS DR16 ELG data, with an enhanced number density (see Table 3) for reducing variances (Avila et al. 2020).

3 METHODOLOGY

In this section we describe the recipe of measuring BAO peak positions from the 2PCFs of LRGs, ELGs, and the corresponding void catalogues, as well as the way for joint constraints using all tracers available. The methods are based on our previous analysis (Zhao et al. 2020), but with various updates. Besides, we briefly introduce the method of extracting cosmological parameters from the BAO measurements.

3.1 BAO reconstruction

Since BAO reconstruction improves the significance of BAO peak in general, for both galaxies and voids (see e.g. Alam et al. 2017; Zhao et al. 2020), we use this technique throughout this work. It is a method introduced by Eisenstein et al. (2007), aiming at sharpening the BAO peak by mitigating nonlinear degradation effects, which are dominated by the bulk flows of large-scale structures and formations of superclusters. This method is then extended (Noh, White & Padmanabhan 2009; Padmanabhan, White & Cohn 2009) and validated using both simulations and observational data (e.g. Seo et al. 2010; Padmanabhan et al. 2012). In brief, the BAO reconstruction algorithm requires an estimation of the dark matter density field based on the distribution of biased tracers, usually in redshift space. This density field is smoothed with a Gaussian kernel, and converted to a displacement field using the Zel’dovich approximation. Finally, displacements are applied to galaxies and random samples to reverse the bulk flow motions. In particular, the shifted random catalogue is used for estimating the galaxy density field after reconstruction (Padmanabhan et al. 2012).

The robustness of the BAO reconstruction technique and the improvement of BAO measurement precision have been verified by numerous studies (e.g. Xu et al. 2012; Burden et al. 2014; Vargas-Magaña et al. 2015; White 2015; Seo et al. 2016). Besides, a number of different reconstruction algorithms are proposed, which aim at better estimations of the displacement fields (e.g. Burden, Percival & Howlett 2015; Hada & Eisenstein 2018; Sarpa et al. 2019; Kitaura et al. 2021; Liu, Yu & Li 2021; Seo et al. 2022). BAO reconstruction is now a standard method for measuring BAO from spectroscopic galaxy data (e.g. Alam et al. 2017; Bautista et al. 2021; Raichoor et al. 2021).

The BOSS DR12 galaxies and the corresponding MD-Patchy mocks are reconstructed using the algorithm described in Padmanabhan et al. (2012). In particular, the reconstructions are run with the full redshift range of 0.2 < z < 0.75, hence the ‘LRG (a)’ and ‘LRG (b)’ samples are reconstructed simultaneously (Alam et al. 2017). Meanwhile, the iterative reconstruction method⁹ (Burden et al. 2015; Bautista et al. 2018) is used for the ‘LRG (c)’ and ‘ELG’ samples, as well as the corresponding EZmock catalogues (Bautista et al. 2021; Raichoor et al. 2021). Note that the CMASS and eBOSS DR16 LRGs are reconstructed as a whole sample, for both the observational data and each mock realization. But the ‘LRG (c)’ and ‘ELG’ samples are reconstructed independently, even though they share the same survey volume to a large extent. For the N-body simulation galaxy catalogues, we rely on the revolver¹⁰ (Nadathur et al. 2019a) implementation of the iterative algorithm for BAO reconstruction.

The input parameters used for the reconstructions of all our samples are detailed in Table 4. In particular, the NGC and SGC samples are reconstructed separately – as they are too far away from each other for cross-correlations to matter – but using the same set of parameters. It is worth noting that small variances of the cosmological parameters and smoothing length do not significantly affect BAO measurements (Vargas-Magaña et al. 2015; Sherwin & White 2019; Carter et al. 2020). Moreover, though the reconstruction algorithms used for different samples in this work are not identical, there are no obvious biases on BAO measurements (see Vargas-Magaña et al. 2018; Bautista et al. 2021; Raichoor et al. 2021, and Section 4). Actually, our reconstructed galaxy samples are the same as those used in eBOSS21.

Table 4.

BAO reconstruction parameters for different samples. Here, f denotes the growth factor for correcting redshift space distortions, b denotes the galaxy bias parameter for estimating the matter density field, N_grid indicates the grid size for sampling the density field, and Σ_r is the Gaussian smoothing length.

Sample	f	b	Ngrid	Σr (⁠|$h^{-1}\, {\rm Mpc}$|⁠)
LRG (a) + (b)	0.757	1.85	5123	15
LRG (c)	0.815	2.3	5123	15
ELG	0.820	1.4	5123	15
MD-Patchy	0.757	2.2	5123	15
EZmock for LRG (c)	0.815	2.3	5123	15
EZmock for ELG	0.820	1.4	5123	15
BigMD (a)	0.716	2.4	5123	10
BigMD (b)	0.757	2.5	5123	10
BigMD (c)	0.797	2.7	5123	10
OuterRim	0.823	1.8	5123	10

3.2 Void catalogue creation

We construct void catalogues upon the reconstructed galaxy samples for both observational data and approximate mocks. They are further selected by the void radii to maximize the BAO significances. We do not generate void catalogues using the N-body simulations, as the void clustering depends on the density and completeness of a galaxy sample, especially for the radial selection function (Forero-Sánchez et al. 2021). In this case, to reproduce the void 2PCFs of the observational data using catalogues from N-body simulations, one has to apply realistic geometric and radial selection effects to the simulated galaxy catalogues. Consequently, the statistical variances of the catalogues generally become too large for precise BAO damping measurements. Thus, we do not constrain the nonlinear damping of void BAO using N-body simulations in this work, but leave it to a future study.

3.2.1 Void finding

We rely on dive for identifying voids from a galaxy sample. It is a fast and parameter-free void finder, based on the cgal¹¹ (The CGAL Project 2021) implementation of the Delaunay triangulation technique, which connects galaxies in comoving space to construct tetrahedra, such that their circumscribed spheres do not contain any galaxy inside. These circumspheres are used as our void tracers, which are essentially the largest spheres in the survey volume that are empty of galaxies. It is worth noting that a high overlapping fraction is found for this type of voids (Zhao et al. 2016), though the associated tetrahedra do not overlap, which turns out to be vital for detecting BAO from voids, as the BAO peak is not prominent for the clustering of disjoint voids (Kitaura et al. 2016a).

Similar to the case of BAO reconstruction, the void finder is run on the full BOSS volume, so voids for the ‘LRG (a)’ and ‘LRG (b)’ samples are generated at once. Besides, the ‘LRG (c)’ void samples are constructed without distinguishing between CMASS and eBOSS LRGs, but independently of the ‘ELG’ voids. Void catalogues generated by dive are then trimmed to fit the survey volume: spheres are kept only if their centres are inside the footprints and redshift ranges, but outside veto masks, of the corresponding samples. For the ‘LRG (c)’ sample, we apply the union of the CMASS and eBOSS LRG footprints. There are however some (<3000) very large spheres with radii larger than |$100\, h^{-1}\, {\rm Mpc}$|⁠, which reside inside the CMASS-only regions, and mostly at z ≳ 0.85, where there are only few galaxies. These objects are removed later during the void radius selection process (see Section 3.2.3), though their impacts on the clustering measurements are negligible on the BAO scale.

The redshift and radius distributions of the resulting void catalogue for each galaxy sample are shown in Fig. 3. It can be seen that the majority of the void radii are between 10 and 20 |$h^{-1}\, {\rm Mpc}$|⁠, which is close the to mean separation of galaxies, but smaller than that of many other void definitions (e.g. Mao et al. 2017; Aubert et al. 2020; Hawken et al. 2020). The is because we define voids purely geometrically, without removing any galaxies, and no merging process is performed even if voids overlap to a large extent. The redshift dependence of void radius distributions is mostly due to the inhomogeneity of the radial number density of galaxies (Forero-Sánchez et al. 2021). For the ‘LRG (c)’ sample, since the number densities of CMASS and eBOSS LRGs differ significantly at high redshifts (z ≳ 0.7, see Fig. 2), we observe two branches in the heat map shown in the middle panel of Fig. 3. These two void populations are found to be located in the overlapping regions between CMASS and eBOSS LRGs, and the CMASS only area (see Fig. 1), respectively.

Figure 3.

Redshift and radius distributions of void samples generated by dive, for each SDSS sample. Colours in the heat maps indicate the comoving number density of voids, and the darker the colour, the higher the void density. Subplots to the right of the heat maps show the mean number density of voids with different radii, in three redshift bins indicated by the colours. The upper histograms show the radial void number density distributions, in radius bins encoded by different colours. Values of number densities are not shown, as the maps and curves are normalized arbitrarily.

Moreover, our spheres with small radii show similar radial distributions as the corresponding galaxies (see Fig. 2), while the distributions of large voids are in general less redshift-dependent, or anticorrelate with that of galaxies. These effects are consistent with the galaxy density dependences of voids in different sizes, studied by Forero-Sánchez et al. (2021) using simulations. The behaviours can be explained by two distinct void populations – voids-in-clouds and voids-in-voids – which reside in regions with overdensities and underdensities, respectively (e.g. Sheth & van de Weygaert 2004; Zhao et al. 2016). A void selection algorithm is then necessary for mitigating contamination of the BAO signature from voids-in-clouds, and maximizing the statistical power of void BAO (Liang et al. 2016). To this end, we have to construct random catalogues for our void samples first to perform clustering measurements.

3.2.2 Random catalogue construction

The clustering signal is defined as the excess spatial correlation of tracers compared to that of a random distribution. Thus a random sample is required for clustering measurements, which should account for effects that may introduce additional clustering pattern, such as the survey geometry and sample completeness in both angular and radial directions. The angular inhomogeneity of galaxies can be assessed by correlating galaxy densities with potential observational systematics, including but not limited to stellar density, seeing, airmass, and Galactic extinction (e.g. Ross et al. 2012, 2020; Raichoor et al. 2021). For the radial selection function, one can either spline interpolate the observed redshift distribution of data to sample redshifts for the random catalogue, or, randomly shuffle redshifts from the data and assign them to the randoms. Both methods show similar clustering measurements based on studies of the BOSS DR9 mocks, and results from the ‘shuffled’ method are slightly closer to the expectations (Ross et al. 2012). Recent studies have shown that identical redshifts in the data and random catalogues may bias the clustering measurements, especially when the survey area is small, which can be corrected using an additional integral constraint term (see de Mattia & Ruhlmann-Kleider 2019; Tamone et al. 2020; Zhao et al. 2021a).

However, sample variations across the survey volume are not trivial for voids, as the number densities of these indirect tracers do not response straightforwardly to observational effects. Besides, the distributions of voids with different radii have to be taken into consideration. In general, the relationship between the void and galaxy number densities depends on the clustering of galaxies, and it is difficult to translate directly the completeness of galaxies to that of voids (Forero-Sánchez et al. 2021). Thus, we follow the random generation scheme described in Liang et al. (2016), which makes use of the ‘shuffled’ method. In this manner, multiple mock realizations are stacked, and the coordinates are shuffled in both angular and radial directions. Moreover, we shuffle void radii along with their redshifts, in different radius and redshift bins, since the void size distributions are sensitive to redshifts, especially for galaxy samples with dramatic variations of radial selection functions (see Fig. 3).

In practice, for each data sample, we stack 100 realizations of mocks catalogues, and divide the combined catalogue into subsamples with different radius and redshift bins. The bin sizes for radii and redshifts are |$1\, h^{-1}\, {\rm Mpc}$| and 0.05 respectively, and adjacent radius bins are merged if there are not enough objects (<1000). Then, we shuffle the (radius, redshift) pairs together for objects inside each bin. Finally, the shuffled subsamples are combined again, and 20 per cent of the objects are randomly chosen to construct the random catalogue for clustering measurements. In this case, the number of random objects are roughly 20 times the number of voids. This is a typical random-to-data ratio for 2PCF estimators (e.g. Vargas-Magaña et al. 2018). Moreover, the chance of having identical redshifts between the data and random catalogues is largely reduced, and given the fact that the biases caused by pairs with identical redshifts are small for the 2PCF monopole (Zhao et al. 2021a), we do not apply integral constraint corrections in this work.

Due to the large angular variations of galaxy densities for the ‘LRG (c)’ sample, depending on whether the CMASS and eBOSS LRG samples overlap, we apply the ‘shuffled’ random generation procedure to voids in the CMASS-only, eBOSS-only, and the overlapping regions individually, and combine the three random catalogues afterwards. Similarly, random samples for ELG chunks – angular regions in which the plate and fibre assignments are performed independently (Raichoor et al. 2021) – are generated separately as well.

3.2.3 Optimal radius selection

With the void catalogues and the corresponding randoms, we are able to measure BAO from the 2PCFs of voids with different radii, using the estimators described in Section 3.3. In order to maximize the cosmological gain from void BAO, we have to avoid contamination from voids-in-clouds that actually reside in overdensity regions, and keep as many voids-in-voids as possible to reduce the statistical uncertainties. To this end, we perform a grid search of the optimal void radius selection criteria, based on the signal-to-noise ratio (SNR) of the void BAO peak from mocks.

Following Liang et al. (2016), we use a model-independent definition of the BAO peak signal S_BAO, which is essentially the height of the BAO peak with respect to the left and right dips:

$$\begin{eqnarray*} S_{\rm BAO} = \left|\, \xi (s_{\rm peak}) - \frac{ \xi (s_{\rm left}) + \xi (s_{\rm right}) }{ 2} \, \right| , \end{eqnarray*}$$

(6)

where ξ(s) denotes the 2PCF monopole at separation s, and s_left, s_peak, and s_right indicate the bins at the left dip, the peak, and the right dip near the BAO scale (⁠|$\sim 100\, h^{-1}\, {\rm Mpc}$|⁠), respectively. The absolute value is necessary since the BAO ‘peak’ of the galaxy–void cross-correlation function is typically negative. We compute 2PCFs with the separation range of |$s \in (0, 200)\, h^{-1}\, {\rm Mpc}$| in this work, with the bin size of |$5\, h^{-1}\, {\rm Mpc}$|⁠. The separation ranges we choose for s_left, s_peak, and s_right are then (80,90), (100,105), and |$(115,125)\, h^{-1}\, {\rm Mpc}$|⁠, respectively. Note that there are two bins for each dip, and we take the average 2PCF value of these bins for ξ(s_left) and ξ(s_right), as is done in Liang et al. (2016). The noise of the BAO signature is then measured as the standard deviation of S_BAO, evaluated from 100 mock realizations.

It has been shown by Forero-Sánchez et al. (2021) that voids selected by constant radius cuts provide near-optimal BAO SNR for BAO reconstructed samples. Moreover, they are not sensitive to moderate observational systematics that introduces <20 per cent incompleteness. Thus, we explore only constant minimum radii R_{v, min} for void selections in this work. In particular, we examine 26 R_{v, min} values, from 12.5 to |$25\, h^{-1}\, {\rm Mpc}$|⁠, with the step size of |$0.5\, h^{-1}\, {\rm Mpc}$|⁠. The void exclusion effect introduces strong anticorrelations on scales below twice the void radii, for both void–void (ξ_vv) and galaxy–void (ξ_gv) 2PCFs (Liang et al. 2016). To prevent its contamination on the BAO peak, which is at around |$100\, h^{-1}\, {\rm Mpc}$| given our current understanding of the Universe (e.g. eBOSS21), we set a maximum radius of |$R_{\rm v, max} = 40\, h^{-1}\, {\rm Mpc}$| for the selections of all our void samples.

The BAO SNR of ξ_vv and ξ_gv for voids with radii between R_{v, min} and R_{v, max}, measured from 100 approximate mock realizations, are shown in Fig. 4, together with the void radius distributions for different samples. The optimal radius cuts for ξ_vv and ξ_gv are in general consistent, and are slightly larger than the most probable void radii, which are at the peaks of the void radius distributions. Hence, we choose a single radius cut for each sample, which yields near-optimal BAO SNR for both ξ_vv and ξ_gv. It is worth noting that we use directly the radius cuts obtained from the ‘LRG (c)’ and ‘ELG’ samples separately, for their cross-correlations (see Section 3.3). These radius cuts, as well as the corresponding void sample sizes, are listed in Table 5. It can also be seen that with the minimum radius cut, we remove 60–70 per cent voids, and the number of remaining voids is about twice the number of the corresponding galaxies (see Table 1) for each sample.

Figure 4.

SNR of the BAO peak, measured from 100 realizations of mock catalogues, for void–void (ξ_vv), galaxy–void (ξ_gv), and galaxy–galaxy (ξ_gg) correlation functions. In particular, ξ_vv and ξ_gv are measured with different minimum void radii (R_{v, min}), and the same maximum radius of |$R_{\rm v, max} = 40\, h^{-1}\, {\rm Mpc}$|⁠. The vertical dotted lines show the void radius cuts used in this work. Shadowed areas show normalized void radius distributions.

Table 5.

Void radius selection ranges (R_{v, min}, R_{v, max}), and the number of voids with different selection criteria, for each galaxy sample. |$N_{\rm void}^{\rm max}$| indicates the number of voids with radii in the range of (0, R_{v, max}), while N_void indicates the number of voids we use for the BAO analysis in this paper, i.e. with radii R_v ∈ (R_{v, min}, R_{v, max}).

Sample	(Rv, min, Rv, max)	Galactic	|$N_{\rm void}^{\rm raw}$|	Nvoid
	|$(h^{-1}\, {\rm Mpc})$|	cap
LRG (a)	(16, 40)	NGC	2639 522	902 637
		SGC	1047 500	317 737
		Total	3687 022	1220 374
LRG (b)	(16, 40)	NGC	3107 491	995 763
		SGC	1130 310	359 202
		Total	4237 801	1354 965
LRG (c)	(21.5, 40)	NGC	1485 650	551 265
		SGC	697 680	243 530
		Total	2183 330	794 795
ELG	(16, 40)	NGC	370 356	172 208
		SGC	362 127	149 988
		Total	732 483	322 196

Furthermore, the optimal BAO SNR of ξ_gv can be higher than that of the galaxy auto 2PCFs (ξ_gg), which are shown as horizontal dashed lines in Fig. 4. This is consistent with the finding in Forero-Sánchez et al. (2021), and suggests that ξ_gv may yield better BAO constraints than ξ_gg. Thus, the galaxy–void cross-correlation is a promising probe of cosmological parameters, considering the fact that it is also less sensitive to observational systematics than galaxy auto correlations (Forero-Sánchez et al. 2021).

Note that the radius selection is the only degree of freedom for the construction of our void sample. Even if the SNR estimated using equation (6) does not represent faithfully the uncertainty of BAO scale determination, there should be no bias for the BAO measurements from voids, though the statistical errors may not be minimized. In fact, we have checked the BAO measurements from the mean 2PCF of mocks with different R_{v, min} values that are close to our fiducial choice, and no significant bias is found. It is also important to stress that apart from the redshift to distance conversion process, which is unavoidable for 3D clustering analysis, the whole void sample construction pipeline is model independent.

3.3 Correlation function estimator

2PCF quantifies the spatial distribution of objects on different scales. In practice, they are typically estimated through pair counts, which are histograms of tracer pairs binned by pair separations. Here, the tracers can be galaxies, voids, and randoms. In this work, all the pair counts are evaluated in 3D comoving space, using the Fast Correlation Function Calculator¹² [fcfc; Zhao in preparation]. Besides, we consider only the monopole 2PCFs, i.e. 2PCFs in isotropic separation bins.

Let |$\widehat{\rm XY} (s)$| be the weighted number of tracer pairs with separation bin s, where members of the pairs are from catalogues X and Y respectively, the normalized pair count is given by

$$\begin{eqnarray*} {\rm XY}(s) = \widehat{\rm XY}(s) / (n_{_{\rm X}} n_{_{\rm Y}}) . \end{eqnarray*}$$

(7)

Here, |$n_{_{\rm X}}$| and |$n_{_{\rm Y}}$| indicate the total weighted number of objects in X and Y, respectively. Note that the denominator of this equation should be replaced by |$n_{_{\rm X}} (n_{_{\rm X}} - 1)$| for auto pair counts, or the case of X = Y. But when |$n_{_{\rm X}} \gg 1$|⁠, equation (7) holds for auto pair counts as well. We shall use this approximation throughout this work, as our numbers of tracers are sufficiently large (see Tables 1 and 5). In particular, we apply the weight expressed by equation (2) to galaxies, but no weight (or a weight of 1) to voids, for the pair counting and normalization processes.

The 2PCF of a sample given the data (D) and random (R) catalogues can then be computed using the Landy–Szalay (LS) estimator (Landy & Szalay 1993):

$$\begin{eqnarray*} \xi = \frac{ {\rm DD} - 2 {\rm DR} + {\rm RR}}{{\rm RR}} . \end{eqnarray*}$$

(8)

Here, we have dropped the separation s for simplicity. This estimator applies to our void auto correlation functions. But for a reconstructed galaxy sample, due to the existence of the shifted random catalogue (S), which is generated by moving random objects with the displacement field for reversing galaxy motions, the 2PCF estimator is slightly different (Padmanabhan et al. 2012):

$$\begin{eqnarray*} \xi _{\rm gg} = \frac{{\rm DD} - 2 {\rm DS} + {\rm SS}}{{\rm RR}} , \end{eqnarray*}$$

(9)

where D denotes the reconstructed galaxy catalogue, and R indicates the unshifted random sample. This equation is further generalized following the spirit of Szapudi & Szalay (1997), for the estimator of the galaxy–void cross-correlation:

$$\begin{eqnarray*} \xi _{\rm gv} = \frac{{\rm D_g D_v - D_g R_v - S_g D_v + S_g R_v}}{{\rm R_g R_v}} . \end{eqnarray*}$$

(10)

Here, the subscripts ‘g’ and ‘v’ indicate catalogues for galaxies and voids, respectively. Furthermore, cross-correlations between ‘LRG (c)’ and ‘ELG’ tracers are estimated using similar formulae.

Usually, we need to combine the 2PCFs of different data sets, such as those for NGC and SGC. This can be done by weighting 2PCFs with the volumes of the galactic caps (or equivalently areas as the redshift ranges are identical), as they are generally too far away from each other to have cross-correlations on scales of interest. However, we combine different samples at the pair count level in this work, as the cross-correlations can be considered for tracers residing in the same volume, with optionally additional weights applied to the samples. This is crucial for the combination of galaxy and void correlations, as is done in Zhao et al. (2020).

To combine two data sets – indicated by subscripts ‘1’ and ‘2’ respectively, and both contain data and a random catalogue – with constant weights, we need only adding weights to one of them, as one can always renormalize weights applied to the two data sets simultaneously, without altering the clustering. Assume that a constant weight w is to be applied to the second data set, given equation (7), the combined normalized data–data pair count can be expressed as

$$\begin{eqnarray*} {\rm DD_{comb}} = \frac{ n_1^2 \cdot {\rm D_1 D_1} + w n_1 n_2 \cdot {\rm D_1 D_2} + w^2 n_2^2 \cdot {\rm D_2 D_2} }{ (n_1 + w n_2)^2 } , \end{eqnarray*}$$

(11)

where n₁ and n₂ are the number of tracers in the data catalogues of the two data sets, respectively. Note, however, that w is generally not the correct weight for the random catalogue, since the ratios of the number of randoms to that of the data can be different for the two data sets. Thus, we introduce an additional weight for randoms of the second data set, which corrects for the random-to-data ratio:

$$\begin{eqnarray*} w_\alpha = \frac{ n_{\rm r1} / n_1 }{ n_{\rm r2} / n_2 } , \end{eqnarray*}$$

(12)

where n_r1 and n_r2 denote that number of random objects of the two data sets. Consequently, the total weight to be applied to the random catalogue of the second data set is

$$\begin{eqnarray*} w_{\rm r} = w \cdot w_\alpha = \frac{ n_2 n_{\rm r1} }{ n_1 n_{\rm r2} } w , \end{eqnarray*}$$

(13)

and the combined random–random and data–random pair counts are then

$$\begin{eqnarray*} {\rm RR_{comb}} = \frac{ n_{\rm r1}^2 \cdot {\rm R_1 R_1} + w_{\rm r} n_{\rm r1} n_{\rm r2} \cdot {\rm R_1 R_2} + w_{\rm r}^2 n_{\rm r2}^2 \cdot {\rm R_2 R_2} }{ (n_{\rm r1} + w_{\rm r} n_{\rm r2})^2 } , \end{eqnarray*}$$

(14)

$$\begin{eqnarray*} {\rm DR_{comb}} &=& \frac{ n_1 n_{\rm r1} \cdot {\rm D_1 R_1} + w w_{\rm r} n_2 n_{\rm r2} \cdot {\rm D_2 R_2} }{ (n_1 + w n_2) (n_{\rm r1} + w_{\rm r} n_{\rm r2})} \nonumber\\ && + \frac{ w_{\rm r} n_1 n_{\rm r2} \cdot {\rm D_1 R_2} + w n_{\rm r1} n_2 \cdot {\rm R_1 D_2} }{(n_1 + w n_2) (n_{\rm r1} + w_{\rm r} n_{\rm r2})} . \end{eqnarray*}$$

(15)

The formulae of the combined pair counts are general, hence one can simply replace D₁, R₁, D₂, and R₂ by the data and random catalogues for galaxies and voids, including the shifted randoms, to estimate the combined pair counts for terms in equations (8–10). In particular, the NGC and SGC pair counts can be combined with w = 1, in which case the cross pair counts between the two data sets are eliminated. Note that correlation functions we use throughout this work are always with NGC and SGC combined, apart from those for consistency examinations in Appendix  E,.

The resulting galaxy–galaxy, galaxy–void, and void–void correlation functions are shown in Fig. 5. Measurements from the SDSS data and the corresponding approximate mocks are generally consistent on BAO scales, with deviations smaller than 1σ for most cases, though the BAO peak of the approximate mocks are typically not as sharp as that of the N-body simulations. Here, errors of the N-body simulation catalogues are estimated using the jackknife resampling method (e.g. Norberg et al. 2009). The discrepancies between data and mocks at |$s \sim 150\, h^{-1}\, {\rm Mpc}$| are possibly due to uncorrected observational systematics, but they do not bias BAO measurements provided models accounting for the broad-band amplitudes (Ross et al. 2017). The cross-correlations between tracers of the ‘LRG (c)’ and ‘ELG’ samples are shown in Fig. 6. Again, the agreements between data and mocks are reasonably well on the BAO scale, and the BAO peaks of cross-correlations between all different tracers are prominent.

Figure 5.

Galaxy–galaxy, galaxy–void, and void–void 2PCFs for different samples, with NGC and SGCs combined. Dots indicate measurements from the SDSS data, with error bars being the standard deviation of measurements from 1000 realizations of the corresponding approximate mocks (Patchy or EZmock). Red-dashed lines and orange envelopes show the mean and 1σ dispersions of 2PCFs from these mocks. Cyan regions denote jackknife error estimations for N-body simulation galaxy catalogues, including BigMD and OuterRim, which are calibrated with the corresponding data. In particular, the 2PCF of the OuterRim simulation is shifted by α = 0.942 to account for the difference of cosmology models (see Table 2).

Figure 6.

Two-point cross-correlation functions between tracers of the ‘LRG (c)’ and ‘ELG’ samples, including the ones for LRGs and ELGs (⁠|$\xi _{\rm gg}^\times$|⁠), LRG voids and ELG voids (⁠|$\xi _{\rm vv}^\times$|⁠), and LRGs and ELG voids (⁠|$\xi _{\rm gv}^\times$|⁠), as well as LRG voids and ELGs (⁠|$\xi _{\rm vg}^\times$|⁠), with northern and southern galactics combined. Red-dashed lines and orange envelopes show the mean and standard deviations measured from 1000 realizations of EZmock catalogues, and dots indicate measurements from the SDSS data, with error bars being the dispersions of mocks.

We then apply different negative weights to the void samples, and combine the galaxy–galaxy, galaxy–void, and void–void pair counts. This is because a negative weight enhances the BAO significance of the combined galaxy and void sample (Zhao et al. 2020). The results for the weight values of −0.05 and −0.2 are shown in Fig. 7. Note that with a void weight of 0, the combined correlation function is essentially the galaxy–galaxy 2PCF. The contribution of voids to the combined correlation function becomes larger when the absolute value of the weight increase. The combined correlation functions we mention hereafter always refer to the results for the combined galaxy and void samples, unless otherwise stated.

Figure 7.

Correlation functions evaluated by combining galaxy–galaxy, galaxy–void, and void–void pair counts, with constant weight applied to voids. They are equivalent to correlation functions of a joint galaxy and void sample, with the same weight for all voids. The left and right columns show results with the void weight of −0.05 and −0.2, respectively. Dots denotes measurements from the SDSS data, and redshift lines indicate mean results from 1000 realizations of EZmock catalogues. The standard deviations of these mocks are shown as orange envelopes, as well as error bars on the data points.

3.4 BAO fitting

In order to measure the BAO peak positions from 2PCFs, we rely on the template fitting method introduced by Xu et al. (2012), and adapted for voids by Zhao et al. (2020) and Variu et al. (in preparation). Since the clustering measurements of galaxies and voids are strongly correlated (see Appendix  A), it is crucial to use a multitracer BAO fitting scheme that takes into account the cross covariances. To this end, we introduce two multitracer approaches:

• Fit the combined 2PCFs with weights applied to different tracers, as is done in Zhao et al. (2020)

• Fit the stacked 2PCFs of multiple tracers, with all their cross covariances included.

3.4.1 BAO models

The theoretical model we use for BAO fitting is based on a template correlation function ξ_t(s), and (Xu et al. 2012)

$$\begin{eqnarray*} \xi _{\rm model} (s) = B^2 \xi _{\rm t} (\alpha s) + A(s), \end{eqnarray*}$$

(16)

where B is a normalization factor that controls the overall amplitude of the model, and A(s) indicates a polynomial that accounts for the broad-band shape, which consists of three nuisance parameters a₀, a₁, and a₂:

$$\begin{eqnarray*} A(s) = a_0 \, s^{-2} + a_1 \, s^{-1} + a_2 . \end{eqnarray*}$$

(17)

It has been shown that this polynomial term yields unbiased BAO measurements (e.g. Xu et al. 2012; Vargas-Magaña et al. 2014). Lastly, α is the BAO dilation parameter, which quantifies the horizontal shift of the model curve, and is essentially the measurement of the BAO peak position. Since α represents the relative difference between the model and template, it can be converted to a ratio of distance scales (Xu et al. 2012), i.e.

$$\begin{eqnarray*} \alpha = \frac{ D_{_{\rm V}} (z) / r_{\rm d} }{ D_{\rm _V, fid} / r_{\rm d, fid}} , \end{eqnarray*}$$

(18)

where r_d is the sound horizon at the drag epoch, and |$D_{_{\rm V}}$| indicates the volume-averaged angular diameter distance (Eisenstein et al. 2005):

$$\begin{eqnarray*} D_{_{\rm V}} (z) = \left[z D_{_{\rm M}}^2 (z) D_{_{\rm H}}(z) \right]^{1/3} = \left[ D_{_{\rm M}}^2 (z) \frac{ {\rm c} z}{H(z)} \right]^{1/3}, \end{eqnarray*}$$

(19)

with |$D_{_{\rm M}} (z)$| and |$D_{_{\rm H}}(z)$| being the angular diameter distance and the Hubble distance, respectively. Besides, H(z) is the Hubble parameter. The subscript ‘fid’ in equation (18) indicates parameters of the cosmology model used for generating the template ξ_t(s).

This template correlation function is actually the Hankel transform of a template power spectrum P_t(k):

$$\begin{eqnarray*} \xi _{\rm t} (s) = \int \frac{k^2\, {\rm d} k}{2 {\pi }^2} \, P_{\rm t} (k) j_0 (k s) \, {\rm e}^{-k^2 a^2} , \end{eqnarray*}$$

(20)

where j₀ is the zero-order spherical Bessel function of the first kind, and a is a parameter for high-k damping, which reduces numerical instability of the integration. The value of a used in this work is |$2\, h^{-1}\, {\rm Mpc}$|⁠, with which the BAO measurements are unbiased, and more robust against noises of the template power spectrum, compared to the results with smaller a values (Variu et al. in preparation).

The template power spectrum for galaxy auto correlation is generated using the ‘de-wiggled’ model (Xu et al. 2012), i.e.

$$\begin{eqnarray*} P_{\rm t, dw} (k) = \left[ P_{\rm lin} (k) - P_{\rm lin, nw} (k) \right] \, {\rm e}^{-k^2 \Sigma _{\rm nl}^2 / 2} + P_{\rm lin, nw} (k) , \end{eqnarray*}$$

(21)

where P_lin(k) is the linear matter power spectrum, P_{lin, nw} indicates its ‘non-wiggle’ counterpart that is free of BAO wiggles, and Σ_nl denotes the BAO damping parameter. The linear matter power spectra we use for the fits, are generated using the camb¹³ (Lewis, Challinor & Lasenby 2000) software. We use our fiducial cosmology to generate the templates for the SDSS data and approximate mocks, while for the N-body simulation catalogues, their corresponding cosmology models are used (see Table 2). Then, P_{lin, nw}(k) is constructed by smoothing BAO wiggles of the corresponding P_lin, using spline fits.

The ‘de-wiggled’ BAO model, however, does not work well for galaxy–void, void–void, and combined correlation functions (Zhao et al. 2020; Variu et al. in preparation). This is because the void exclusion effect not only introduces high-k oscillatory patterns in the void–void and galaxy–void power spectra, but also alter their broad-band shapes (Chan, Hamaus & Desjacques 2014; Zhao et al. 2016). The polynomial term in equation (17) is not able to describe precisely the complicated shapes. In order to account for these effects, we generalize equation (21) by introducing a term P_{t, nw}(k), that models the broad-band power spectrum shape (see Zhao et al. 2020, for details):

$$\begin{eqnarray*} P_{\rm t} (k) = P_{\rm t, dw} (k) \cdot \frac{ P_{\rm t, nw} (k) }{ P_{\rm lin, nw} (k) } . \end{eqnarray*}$$

(22)

This additional ‘non-wiggle’ term encodes both the broad-band clustering signal and geometric effects that are closely related to the void definition through Delaunay triangulations. Thus, it is difficult to develop an accurate analytical model for this term. To circumvent this problem, we rely on different models for the two multitracer analysis methods:

• Fit the combined 2PCFs using a parabolic model with an additional free parameter, as is done in Zhao et al. (2020).

• Fit the stacked 2PCFs of multiple tracers using a template-based numerical model introduced by Variu et al. (in preparation).

3.4.2 Parabolic model for the combined 2PCFs

It has been shown by Zhao et al. (2020) that the generalized BAO model in equation (22) works for the combined 2PCFs of galaxies and voids. Since the P_{t, nw}(k) term depends on the weight applied to voids, in principle, we have to find a formula that is able to describe well the ‘non-wiggle’ term P_{t, nw}/P_{lin, nw} for all the three correlations: galaxy–galaxy, galaxy–void, and void–void. Fig. 8 shows the measurements of P_t/P_{lin, nw} from mocks for different correlations. Though the BAO wiggles are included, the broad-band shapes are the same as those of the corresponding ‘non-wiggle’ terms (see also Zhao et al. 2020). It can be seen that the low-k parts (⁠|$k \lesssim 0.15\, h\, {\rm Mpc}^{-1}$|⁠) of these measurements are consistent with parabolas. Thus, we use the parabolic approximation of the P_{t, nw}/P_{lin, nw} term introduced by Zhao et al. (2020):

$$\begin{eqnarray*} P_{\rm t, para} (k) \approx P_{\rm t, dw} (k) \cdot (1 + c k^2). \end{eqnarray*}$$

(23)

Here, c is parameter that accounts for the shape of P_{t, nw}/P_{lin, nw} at low k. Even though it does not model well the high-k behaviour of measurements shown in Fig. 8, the model works reasonably well as the Gaussian damping term in equation (20) suppresses the high-k discrepancies to a large extent (see Zhao et al. 2020, and Section 4 for the performance of this model).

Figure 8.

Ratio of the mock or template power spectra, to the linear non-wiggle matter power spectrum. Open symbols indicate the mean results measured from 100 realizations of Patchy or EZmock catalogues, and black lines are obtained using the numerical power spectrum templates without BAO wiggles, which are used for the BAO fitting model in equation (22). Values of the ratios are not shown, as the curves are normalized arbitrarily.

3.4.3 Template-based numerical model for individual 2PCFs

Another way of modelling P_{t, nw}/P_{lin, nw} is to use numerical templates that can be constructed using mocks that are free of BAO wiggles. In particular, it is preferable to avoid model dependence during the generation of these templates, other than that of the matter power spectrum. To this end, we construct mock galaxy catalogues upon Gaussian random fields that are generated with the ‘non-wiggle’ power spectrum P_{lin, nw}(k), followed by survey geometry application, including in particular radial selection functions, for all our samples to be fitted. This is performed with the Cosmological GAussian Mock gEnerator¹⁴ (CosmoGAME; Variu et al. in preparation), which implements a simplified version of the galaxy bias model of EZmock (Chuang et al. 2015a; Zhao et al. 2021a). We then construct the corresponding void catalogues following Section 3.2, and measure both the void–void and galaxy–void power spectra. These power spectra are calibrated to those measured from the Patchy and EZmock catalogues, by tuning the galaxy bias parameters of CosmoGAME. Finally, we generate 3000 realizations of Gaussian galaxy and void catalogues for each of our data sample, thanks to the high efficiency of CosmoGAME. The numerical templates of the P_{t, nw}(k) term in equation (22), denoted by |$\mathcal {P}_{\rm gv}$| and |$\mathcal {P}_{\rm vv}$| for galaxy–void and void–void power spectra, respectively, are then generated by averaging the measurements of all the Gaussian mock realizations. Here, the power spectra are all measured using the powspec code.¹⁵ We do not generate templates for P_{gg, nw}, as the de-wiggled model works already well for galaxy auto correlations.

The resulting |$\mathcal {P}_{\rm gv}$| and |$\mathcal {P}_{\rm vv}$| templates for different samples, normalized by P_{lin, nw}, are shown in Fig. 8, in comparisons with the corresponding power spectra measurements from Patchy and EZmock catalogues. In general, these templates capture the main features of the measurements from the approximate mocks. There is, however, an obvious discrepancy for the ELG void auto power spectrum at small k. This is not a critical problem, as small deviations of the template non-wiggle power spectra do not bias BAO measurements (Variu et al. in preparation). After all, residuals in simple shapes can be accounted for by the polynomial term in equation (17). Note that we do not construct templates specifically for the cross-correlations between the ‘LRG (c)’ and ‘ELG’ samples, as the templates for the individual samples are found to be performing well (see Fig. C2).

In principle, it is possible to measure the template power spectrum for the combined sample of galaxies and voids as well, by combining the Gaussian mocks of galaxies and voids with different void weights. We will investigate the possibility of this approach in a future work.

3.4.4 BAO fitting parameter inference

To obtain constraints on the BAO fitting parameters |$\boldsymbol{p}$|⁠, especially for the α parameter, we assume a multivariate Gaussian likelihood:

$$\begin{eqnarray*} \mathcal {L} \approx {\rm e}^{ - \chi ^2 (\boldsymbol{p}) / 2 } , \end{eqnarray*}$$

(24)

with the chi-squared function given by

$$\begin{eqnarray*} \chi ^2 (\boldsymbol{p}) = \boldsymbol{d}^{\sf T} {\bf C}^{-1} \boldsymbol{d} . \end{eqnarray*}$$

(25)

Here, C indicates the covariance matrix (see Appendix  A), and |$\boldsymbol{d}$| denotes the difference between the data and model vectors:

$$\begin{eqnarray*} \boldsymbol{d} = \boldsymbol{\xi }_{\rm data} - \boldsymbol{\xi }_{\rm model} (\boldsymbol{p}) . \end{eqnarray*}$$

(26)

Furthermore, we use the reduced chi-squared, i.e. the best-fitting chi-squared value normalized by the degree of freedom, for examining the goodness of the fits.

For the SDSS data and approximate mocks, the covariance matrices are estimated using the sample variances of mocks, i.e.

$$\begin{eqnarray*} {\bf C}_{\rm s} = \frac{1}{N_{\rm m} - 1} {\bf M}^{\sf T} {\bf M} , \end{eqnarray*}$$

(27)

where N_m denotes the number of mock realizations, and M indicate differences between the 2PCFs of individual mocks and their mean. Elements of M are given by

$$\begin{eqnarray*} {\bf M}_{i j} = \xi _i (s_j) - \bar{\xi } (s_j) , \quad i = 1, 2, \dots , N_{\rm m}, \quad s_j \in (s_{\rm min}, s_{\rm max}) , \end{eqnarray*}$$

(28)

where ξ_i indicates the 2PCF of the i-th mock, |$\bar{\xi }$| is the mean 2PCF of all mocks, and (s_min, s_max) denotes the separation range used for the fit. For N-body simulation catalogues, the covariance matrices are estimated using the jackknife resampling method. For N_sub subvolumes, with one being removed from the sample at each time, the covariance matrix can be estimated by (Norberg et al. 2009)

$$\begin{eqnarray*} {\bf C}_{\rm jk} = \frac{ N_{\rm sub} - 1 }{ N_{\rm sub} } {{\bf M}^{\prime }}^{\sf T} {\bf M}^{\prime } , \end{eqnarray*}$$

(29)

with M′ computed from the 2PCFs of all N_sub subsamples, following equation (28).

Denoting the number of separation bins as N_bin, the unbiased inverse covariance matrix in equation (25) is (Hartlap, Simon & Schneider 2007)

$$\begin{eqnarray*} {\bf C}^{-1} = (1 - D) \, {\bf C}_{\rm s}^{-1}, \end{eqnarray*}$$

(30)

with the correction factor

$$\begin{eqnarray*} D = \frac{N_{\rm bin} + 1}{ N_{\rm m} - 1 } . \end{eqnarray*}$$

(31)

In addition, to assess the uncertainty propagated from the covariance matrix, the fitted parameter covariances are rescaled by (Percival et al. 2014)

$$\begin{eqnarray*} m_1 = \frac{ 1 + B_0 (N_{\rm bin} - N_{\rm par}) }{ 1 + A_0 + B_0 (N_{\rm par} + 1)} , \end{eqnarray*}$$

(32)

and variances of the best-fitting parameters obtained from mocks used for C_s estimation should be rescaled by

$$\begin{eqnarray*} m_2 = (1 - D)^{-1} m_1 . \end{eqnarray*}$$

(33)

Here, N_par is the number of free parameters for the fits, and the factors are

$$\begin{eqnarray*} A_0 = \frac{2}{ (N_{\rm m} - N_{\rm bin} - 1)(N_{\rm m} - N_{\rm bin} - 4) } , \end{eqnarray*}$$

(34)

$$\begin{eqnarray*} B_0 = \frac{N_{\rm m} - N_{\rm bin} - 2}{ (N_{\rm m} - N_{\rm bin} - 1)(N_{\rm m} - N_{\rm bin} - 4) } . \end{eqnarray*}$$

(35)

In practice, we use the BAO Fitter for muLtI-Tracers¹⁶ (BAOflit) to constrain BAO fitting parameters. It encodes all the BAO models described in Section 3.4.1, and is able to fit a single α parameter for 2PCFs of multitracers, with their cross covariances taken into account. The nuisance parameters in equation (17) are computed analytically using the least squares method; while for the rest of the parameters, BAOflit relies on the MultiNest¹⁷ tool (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2019) for accurate and efficient Monte–Carlo (MC) Bayesian posterior sampling. In fact, including the nuisance parameters in the MC sampler does not change the fitting results, but slows down the sampling process significantly due to the enlargement of the parameter space dimension (see also Zhao et al. 2020). The best-fitting (median) value and 1σ confidence limits of the parameters – defined as the 16th (σ₋) and 84th (σ₊) percentiles of the cumulative posterior distribution respectively – are computed using the pymultinest package¹⁸ (Buchner et al. 2014).

Moreover, BAOflit maximizes the numerical stability of matrix manipulations, and speeds up the chi-squared function evaluations, including the least squares calculations, using the QR decomposition technique (see Appendix  B). This optimization can be used for all fitting applications that rely on covariance matrices from mocks or jackknife resampling, especially for those with nuisance parameters that can be calculated using the least squares method.

3.5 Cosmological parameter constraints

We consider in this work the standard flat-ΛCDM model, as well as its two one-parameter extensions, oΛCDM and flat-wCDM, which add the Ω_k and w parameters, respectively. In particular, the Hubble parameter H(z) in these models are given by

$$\begin{eqnarray*} H^2 (z) = H_0^2 \left[ \Omega _{\rm m} (1+z)^3 + \Omega _{\rm k} (1+z)^2 + \Omega _{\rm de} (1+z)^{3(1+w)} \right] , \end{eqnarray*}$$

(36)

where |$H_0 = H(z = 0) \equiv 100\, h$|⁠, and Ω_m + Ω_k + Ω_de = 1. Here, Ω_k and Ω_de denote the spatial curvature density and dark energy density at z = 0 respectively, and w is the dark energy equation of state. When Ω_k = 0 and w = −1, the two extended models reduce to the flat-ΛCDM model.

We compare the parameter constraints with |$D_{_{\rm V}} /r_{\rm d}$| measured in this work via the multitracer BAO analysis, and the corresponding BAO-only results in eBOSS21. Note that in the latter study, the anisotropic BAO measurements expressed by |$D_{_{\rm M}} /r_{\rm d}$| and |$D_{_{\rm H}} /r_{\rm d}$| are used whenever possible, to maximize the constraining power. To break the degeneracies of cosmological parameters measured from BAO, we introduce also the primordial deuterium abundance (Cooke, Pettini & Steidel 2018) for big bang Nucleosynthesis (BBN) calculations, the Planck CMB temperature and polarization data (Planck Collaboration 2020), and the Pantheon type Ia supernovae (SNe Ia) sample (Scolnic et al. 2018), and investigate different combinations of these probes.

In practice, we use a variant¹⁹ of the cosmomC package²⁰ (Lewis & Bridle 2002) for cosmological fittings and model evaluations with different probes throughout this work. For all cases, we run four parallel MC Markov chains, with the first 30 per cent samples removed as the burn in period. The Gelman and Rubin convergence statistics is chosen to be R − 1 < 10⁻³, to have accurate confidence limits for comparing results from different BAO measurements. The posterior distribution plots and marginalized statistics are generated using the getdist tool²¹ (Lewis 2019).

4 ROBUSTNESS AND SYSTEMATIC ERROR ANALYSIS USING MOCK CATALOGUES

In order to investigate the robustness of our analysis, and assess systematic uncertainties, we perform extensive tests of the multitracer BAO constraint methodology, using the large sets of approximate mocks. As stated before, we apply two multitracer BAO fitting methods, including a simple fit to the combined correlation functions evaluated using the pair counts for different tracers, as well as a joint fit with the same α for all individual measurements. We examine the performances of these two approaches using mocks, and choose the more reliable one for the SDSS data analysis. To this end, we compare results with different fitting schemes, such as different fitting ranges and priors of parameters, to estimate potential biases of the measurements.

The fiducial BAO fitting schemes for different samples, as listed in Table 6, are chosen based on these studies. When different clustering statistics are fitted simultaneously with the same α parameter, their own fitting methods are used for the corresponding segments of the data vector. The covariance matrices we use for the studies are all estimated from the same sets of mocks (see Appendix  A). For simplicity, the fitted error on α from the mocks are often quoted as a single value: σ_α = (σ_{α, +} − σ_{α, −})/2. This is because the posterior distributions of α from the mocks are typically highly symmetric. In this case, the median value of the posterior can be estimated as α_{fit, med} ≈ (σ_{α, +} + σ_{α, −})/2.

Table 6.

Fiducial BAO fitting scheme, including the theoretical model, fitting range, and parameter priors, for different samples in this work. For the parameter priors, single numbers indicate fixed values for the fits, while ranges of numbers indicate the extents of flat priors. Σ_{nl, prior} values in parenthesis denote best-fitting values from N-body simulations, while the rest of the figures are all from approximate mocks. Optimal void weights used for the combined correlation functions are also shown, as the superscripts of ξ_comb. In particular, for the cross-correlations between ‘LRG (c)’ and ‘ELG’, the optimal weights for LRG voids and ELG voids are both −0.1. Fitting schemes for the multitracer SDSS data analysis are highlighted in bold.

Sample	Clustering	Model	Fitting range (⁠|$h^{-1}\, {\rm Mpc}$|⁠)	Bprior	Σnl, prior (⁠|$h^{-1}\, {\rm Mpc}$|⁠)	cprior (⁠|$h^{-2}\, {\rm Mpc}^2$|⁠)
LRG (a)	ξgg	(21)	(50, 160)	(1.6, 1.8)	4.9 (2.5)	–
	ξgv	(22) + |$\mathcal {P}_{\rm gv}^{_{\rm LRG (a)}}$|	(60, 160)	(6.2,7.0)	6.1	–
	ξvv	(22) + |$\mathcal {P}_{\rm vv}^{_{\rm LRG (a)}}$|	(70, 160)	(3.3,3.9)	6.9	–
	|$\boldsymbol{\xi }_{\bf comb}^{\boldsymbol{w}{\bf =-0.05}}$|	(23)	(60, 150)	(2.6, 3.0)	(0, 10)	160
LRG (b)	ξgg	(21)	(50, 160)	(1.6,1.8)	4.7 (4.0)	–
	ξgv	(22) + |$\mathcal {P}_{\rm gv}^{_{\rm LRG (b)}}$|	(60, 160)	(6.1,7.4)	6.3	–
	ξvv	(22) + |$\mathcal {P}_{\rm vv}^{_{\rm LRG (b)}}$|	(70, 160)	(3.3,4.3)	7.1	–
	|$\boldsymbol{\xi }_{\bf comb}^{\boldsymbol{w}{\bf =-0.05}}$|	(23)	(60, 150)	(2.5, 2.9)	(0, 9.4)	150
LRG (c)	ξgg	(21)	(50, 160)	(1.5,1.7)	4.8 (4.0)	–
	ξgv	(22) + |$\mathcal {P}_{\rm gv}^{_{\rm LRG (c)}}$|	(60, 160)	(9.5,11)	9.5	–
	ξvv	(22) + |$\mathcal {P}_{\rm vv}^{_{\rm LRG (c)}}$|	(70, 160)	(7.9,8.5)	13	–
	|$\boldsymbol{\xi }_{\bf comb}^{\boldsymbol{w}{\bf =-0.1}}$|	(23)	(60, 150)	(4.1, 4.8)	(0, 11)	13
ELG	ξgg	(21)	(50, 160)	(0.86,0.97)	3.7 (1.5)	–
	ξgv	(22) + |$\mathcal {P}_{\rm gv}^{_{\rm ELG}}$|	(60, 160)	(4.1,5.6)	6.6	–
	ξvv	(22) + |$\mathcal {P}_{\rm vv}^{_{\rm ELG}}$|	(70, 160)	(3.2,3.9)	10	–
	|$\boldsymbol{\xi }_{\bf comb}^{\boldsymbol{w}{\bf =-0.1}}$|	(23)	(60, 150)	(1.5, 2.0)	(0, 13)	70
LRG (c) × ELG	|$\xi _{\rm gg}^\times$|	(21)	(50, 160)	(1.1,1.4)	6.0	–
	|$\xi _{\rm gv}^\times$|	(22) + |$\mathcal {P}_{\rm gv}^{_{\rm ELG}}$|	(60, 160)	(5.5,7.0)	7.9	–
	|$\xi _{\rm vg}^\times$|	(22) + |$\mathcal {P}_{\rm gv}^{_{\rm LRG (c)}}$|	(60, 160)	(5.4,6.3)	7.7	–
	|$\xi _{\rm vv}^\times$|	(22) + |$\mathcal {P}_{\rm vv}^{_{\rm LRG (c)}}$|	(70, 160)	(3.0,3.6)	7.9	–
	|$\boldsymbol{\xi }_{\bf comb}^{\boldsymbol{\times }, \boldsymbol{w}{\bf =\lbrace -0.1,-0.1\rbrace }}$|	(23)	(60, 150)	(1.6, 2.1)	(0, 14)	320

4.1 Optimal void weight

The negative BAO ‘peak’ of the galaxy–void cross-correlation (see Figs 5 and 6) suggests that galaxies and voids are actually anticorrelated on the BAO scale. Therefore, the BAO signal will be diluted if combining galaxies and voids by merging the catalogues directly. Instead, a negative weight can be applied to one of the tracers, to maximize the BAO significance of the combined sample (Zhao et al. 2020). Actually, weights are generally needed when combining multiple tracers with different biases, to optimize the clustering measurements (see e.g. Percival, Verde & Peacock 2004). To this end, we measure the combined correlation functions with different void weights for pair count combinations (see Section 3.3) and explore the optimal weight that minimizes the fitted error of α. In particular, we focus on negative void weights, and examine weights from 0 down to −0.2, with a step size of 0.01. This range has already been shown to be sufficiently wide for the optimal w value determination with BOSS LRGs (cf. Zhao et al. 2020; Forero-Sánchez et al. 2021, for the BOSS results with larger w ranges).

We then fit the mean combined correlation functions of all 1000 mock realizations with the fitting range of |$(60,150)\, h^{-1}\, {\rm Mpc}$|⁠, which is chosen based on Appendix  C. Here, the covariances are rescaled by 1/1000 to allow for tight constraints on all fitting parameters, which are used later for the prior investigations in Appendix  D. The resulting median values of α from the posterior distributions, and the corresponding fitted errors, are shown in Fig. 9. The optimal void weights we choose for the ‘LRG (a)’ and ‘LRG (b)’ samples are both −0.05, which is consistent with the results in Zhao et al. (2020), while for ‘LRG (c)’ and ‘ELG’ the chosen values are −0.1. With these weights, the best-fitting α values of the combined correlation functions are consistent with those of the corresponding galaxy auto correlations (i.e. the combined 2PCF with w = 0), while the fitted errors are significantly reduced. The difference between the measured α and the expected value is the largest for the ‘ELG’ sample, possibly because of the observational systematics applied to the mocks (Raichoor et al. 2021; Zhao et al. 2021a).

Figure 9.

The best-fitting (median) value and 1σ error of α, obtained from fits to the mean combined correlation functions of all mocks, evaluated with different void weights, for each sample. The best-fitting values are shown as relative differences to α_exp, which indicates the expected α value given the cosmology models of the mocks. N_m = 1000 denotes the total number of mock realizations. Vertical-dotted lines indicate the void weights for different samples used in this work.

For the cross-correlations between the ‘LRG (c)’ and ‘ELG’ sample, in principle the optimal weights of the LRG voids (⁠|$w_{_{\rm LRG}}$|⁠) and ELG voids (⁠|$w_{_{\rm ELG}}$|⁠) can be different. For simplicity, we vary in each time the weight for one void sample only, and fix the other void weight to −0.1, which is the optimal value for both individual samples. The fitting results are shown in Fig. 10, where the fitted error is minimized when both void weights are −0.1. The bias of α is smaller than that of the ‘ELG’ sample. This can be explained by the fact that the cross-correlation suffers less observational systematic effects (e.g. Wang et al. 2020; Zhao et al. 2021a). The combined correlation functions discussed hereafter are always the ones with the optimal void weights shown in Table 6, unless otherwise stated.

Figure 10.

Same as Fig. 9, but for the cross-correlation between the ‘LRG (c)’ and ‘ELG’ samples. The left-hand panel shows results with a fixed ELG void weight of −0.1, when varying the weight for LRG voids. In contrast, the right-hand panel is for different ELG void weights, with a fixed LEG void weight of −0.1.

4.2 Fitting ranges

The choice of the fitting range depends on scales that are useful for extracting the BAO peak position, and how well the correlation functions on different scales are reproduced by the templates and broad-band terms of our theoretical models. To understand how they affect our BAO measurements, we run BAO fits with different fitting ranges, to the mean 2PCFs computed using 1000 realizations of mocks, for all the samples. Again, the covariance matrices are divided by 1000. The BAO measurements with different fitting ranges for individual clustering measurements, as well as the combined 2PCFs, are presented in Appendix  C. Actually, our fiducial fitting ranges are chosen based on these results.

We then fit multiple correlation functions simultaneously with the same α, using the corresponding fiducial fitting ranges of individual clustering statistics. In particular, since the fitting results for void–void correlations are sometimes sensitive to the fitting ranges (see Appendix  C), including the void auto correlations (ξ_vv), and the cross-correlation between LRG voids and ELG voids (⁠|$\xi _{\rm vv}^\times$|⁠), we present two sets of joint fits that differ in whether the void–void correlation is included. The fitting results for different samples are shown in Table 7, together with measurements from the combined correlation functions. Here, the fitted errors on α rescaled by |$\sqrt{1000}$|⁠, for approximating the precision of BAO measurements with the observational data.

Table 7.

Fitting results for the mean 2PCFs of 1000 approximate mock realizations, with the covariances evaluated from the same set of mocks, but rescaled by 1/1000. When a data vector consists of multiple clustering statistics, the 2PCFs share the same α parameter during the fit. α_exp denotes the expected α value for each sample, given the cosmology models for constructing the mocks, as well as the fiducial cosmology used for coordinate conversions and BAO fits (see Table 2). α_fit indicates results with the fiducial fitting ranges – [s_smin, s_max], listed in Table 6 – with the fitted errors rescaled by |$\sqrt{1000}$|⁠. N_{bin, fid} indicates the length of the data vector with the fiducial fitting range. |$\Delta _{s_{\rm min}, s_{\rm max}} \alpha _{\rm fit, med}$| denotes the differences of the fitted median values of α, when varying s_min or s_max by |$5\, h^{-1}\, {\rm Mpc}$|⁠.

Sample	αexp	Data vector	Nbin, fid	αfit	αfit, med − αexp	|$\max | \Delta _{s_{\rm min}, s_{\rm max}} \alpha _{\rm fit, med} |$|
MD-Patchy LRG (a)	0.9993	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	22	|$0.9999_{-0.0127}^{+0.0124}$|	0.0005	0.0004
		|$\lbrace \boldsymbol{\xi }_{\rm gv}\rbrace$|	20	|$0.9995_{-0.0121}^{+0.0119}$|	0.0002	0.0002
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	42	|$1.0003_{-0.0111}^{+0.0114}$|	0.0010	0.0003
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	60	|$1.0007_{-0.0113}^{+0.0111}$|	0.0014	0.0007
		|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	18	|$0.9994_{-0.0114}^{+0.0112}$|	0.0001	0.0001
MD-Patchy LRG (b)	0.9996	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	22	|$1.0017_{-0.0118}^{+0.0111}$|	0.0020	0.0002
		|$\lbrace \boldsymbol{\xi }_{\rm gv}\rbrace$|	20	|$1.0006_{-0.0113}^{+0.0118}$|	0.0009	0.0000
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	42	|$1.0017_{-0.0105}^{+0.0103}$|	0.0021	0.0002
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	60	|$1.0014_{-0.0106}^{+0.0106}$|	0.0018	0.0005
		|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	18	|$1.0007_{-0.0105}^{+0.0105}$|	0.0011	0.0001
EZmock LRG (c)	1.0000	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	22	|$1.0012_{-0.0151}^{+0.0154}$|	0.0012	0.0004
		|$\lbrace \boldsymbol{\xi }_{\rm gv}\rbrace$|	20	|$0.9988_{-0.0169}^{+0.0168}$|	−0.0012	0.0004
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	42	|$1.0011_{-0.0143}^{+0.0145}$|	0.0011	0.0005
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	60	|$0.9993_{-0.0135}^{+0.0138}$|	−0.0007	0.0022
		|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	18	|$1.0019_{-0.0145}^{+0.0137}$|	0.0019	0.0016
EZmock ELG	1.0003	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	22	|$1.0059_{-0.0409}^{+0.0395}$|	0.0056	0.0005
		|$\lbrace \boldsymbol{\xi }_{\rm gv}\rbrace$|	20	|$1.0021_{-0.0428}^{+0.0419}$|	0.0018	0.0013
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	42	|$1.0057_{-0.0377}^{+0.0374}$|	0.0054	0.0019
		|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	60	|$1.0065_{-0.0374}^{+0.0384}$|	0.0062	0.0008
		|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	18	|$1.0068_{-0.0383}^{+0.0349}$|	0.0066	0.0020
EZmock LRG (c) × ELG	1.0001	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times \rbrace$|	22	|$1.0023_{-0.0375}^{+0.0385}$|	0.0021	0.0001
		|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times \rbrace$|	62	|$1.0022_{-0.0371}^{+0.0369}$|	0.0021	0.0013
		|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times ,\boldsymbol{\xi }_{\rm vv}^\times \rbrace$|	80	|$1.0017_{-0.0389}^{+0.0405}$|	0.0015	0.0039
		|$\lbrace \boldsymbol{\xi }_{\rm comb}^\times \rbrace$|	18	|$1.0006_{-0.0377}^{+0.0376}$|	0.0005	0.0005

In general, results from the galaxy–void cross-correlations (ξ_gv) and galaxy–galaxy auto correlations (ξ_gg) are similar, in terms of both the median values and fitted errors. Including galaxy–void cross-correlations in the data vectors reduces the BAO measurement uncertainties, compared to the results from galaxies alone. However, it does not necessarily help when further including the void–void correlation. Results from the combined correlation functions are generally consistent with the ones from joint fits to |$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$| and |$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times ,\boldsymbol{\xi }_{\rm vv}^\times \rbrace$|⁠. For all cases, no obvious biases of α are observed, when taking the differences between the fitted and expected α values of the mocks, given the uncertainties of α measurements. It reveals that with the BAO fitting methods and reconstruction algorithms used in this work, we are able to recover the ‘true’ cosmology of the mocks through both multitracer BAO fitting methods, despite of the observational systematics encoded in the catalogues.

In order to assess the potential biases of α due to the fitting range choices, we further vary either the lower or upper bound of the fitting ranges by one separation bin of the 2PCFs (⁠|$5\, h^{-1}\, {\rm Mpc}$|⁠), resulting in four additional fitting range settings. The systematic uncertainties due to the fitting ranges are then estimated as the maximum difference between the resulting α of these four cases, and that of the fiducial fitting range. The biases measured in this way are not significant as well, as is shown in Table 7.

4.3 Parameter priors

We introduce priors of the BAO fitting parameters based on the results from the mean 2PCFs of mocks. The priors are then verified by comparing the fitted errors and the dispersion of α fitted from individual mocks. In particular, we measure the distribution of the pull quantity of α (see e.g. Xu et al. 2012; Bautista et al. 2021),

$$\begin{eqnarray*} g(\alpha) = \frac{ \alpha - \bar{\alpha } }{\sigma _{\alpha }} , \end{eqnarray*}$$

(37)

and compare it to a standard normal distribution. Here, α and σ_α denote the fitted median and 1σ error of α respectively, for each mock realization. |$\bar{\alpha }$| indicates the averaged best-fitting α over all realizations.

Throughout this work, we use a flat prior of α, in the range of (0.8, 1.2). The priors of the rest of the parameters are listed in Table 6, based on investigations in Appendix  D. In summary, we use always flat priors on the B parameter in equation (16), with ranges chosen based on the pull quantities obtained from fits to individual mocks. Since the BAO damping of approximate mocks are generally overestimated compared to the results from N-body simulations (see Table: 6), flat Σ_nl priors are used for the combined correlation functions, with a lower limit of 0, and upper limits estimated using approximate mocks. Even though this type of Σ_nl prior works well for the individual clustering measurements, as well as the combined correlation function, the fitted errors are underestimated for joint fits to stacked 2PCFs (see Table D3). Therefore, we follow the traditional BAO fitting scheme (e.g. Xu et al. 2012; Alam et al. 2017; Raichoor et al. 2021), and fix the Σ_nl parameter when fitting multiple 2PCFs simultaneously. Lastly, for the combined 2PCFs, the c parameter in equation (23) is always fixed, with the value measured from the mean 2PCFs of mocks. Given these priors, the difference between the dispersion of best-fitting α and the median fitted error are mostly less than 5 per cent (see Tables D2 and D3).

With the fiducial choices of priors, the distributions of g(α) for different samples are shown in Fig. 11. They are generally in good agreements with a standard normal distribution. Besides, the distributions of the fitted median values of α are shown in Fig. 12. The histograms are well consistent with the posterior distributions of α from fits to the mean 2PCFs of mocks, i.e. the ones for the results in Table 7. Here, the widths of the posterior distributions are rescaled by |$\sqrt{1000}$|⁠, with respect to the median values. These results validate that our measurements with the fiducial priors are robust, including both best-fitting values and error estimates.

Figure 11.

Distributions of the pull statistic defined in equation (37), based on fits to individual mock realizations, with parameter priors given in Table 6. Black lines shown standard normal distributions, and histograms show results with different fitting data vectors. ‘gg’, ‘gg + gv’, and ‘gg + gv + vv’ denote results from fits to |$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|⁠, |$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|⁠, and |$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$| respectively, for the individual samples. While for the cross-correlations between ‘LRG (c)’ and ‘ELG’, they are for |$\lbrace \boldsymbol{\xi }_{\rm gg}^\times \rbrace$|⁠, |$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times \rbrace$|⁠, and |$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times ,\boldsymbol{\xi }_{\rm vv}^\times \rbrace$|⁠, respectively. ‘comb’ indicates results for the combined correlation functions, with optimal void weights.

Figure 12.

Distributions of the best-fitting (median) values of α. Scatter plots show comparisons of measurements from individual mocks, with different data vectors (see Fig. 11 for the interpretations of subscripts), and the corresponding histograms are shown in red. Black solid curves show posterior distributions of α, obtained from fits to the mean 2PCFs of mocks, with the widths rescaled by |$\sqrt{1000}$|⁠. Purple plus symbols indicate median values drawn from the black curves. Circles and dashed lines in blue denote results from SDSS data.

To assess potential systematic uncertainties due to the choices of parameter priors, we vary our fiducial priors and perform fits to individual mock realizations. In practice, we exchange the prior types used for fitting multiple 2PCFs and the combined correlation function. That is to say, we use flat priors of Σ_nl for joint fits to multiple 2PCFs, and fixed Σ_nl values for the combined correlation functions, with ranges and values taken from Tables D2 and D3. In addition, to examine possible biases due to the c parameter for fitting combined 2PCFs, we vary the c parameter values in Table 6, by 10 per cent for the BOSS tracers and 30 per cent for the eBOSS tracers, respectively, while keeping the fiducial priors of B and Σ_nl. These fractions are roughly the relative errors on c, from fits to the mean 2PCFs of mocks. The biases of the median value of best-fitting α from all mocks, compared to the results with the fiducial priors, are shown in Table 8. The results suggest that our BAO measurements are not sensitive to the prior types and ranges, as the biases are small, especially for the multitracer approach with combined 2PCFs.

Table 8.

Fitting results of individual mock realizations with different parameter priors. 〈α_fid〉 denotes the median value of best-fitting α for all mocks, with the fiducial priors in Table 6. |$\Delta _{B,\Sigma _{\rm nl}} \langle \alpha \rangle$| indicates the bias of 〈α〉 when using a different set of priors for B and Σ_nl, taken from Tables D1 and D2. Δ_c〈α〉 denotes the differences of 〈α〉, when increasing or reducing the fiducial c parameter value, by 10 per cent for the BOSS tracers and 30 per cent for the eBOSS tracers.

Sample	Data vector	〈αfid〉	|$|\Delta _{B,\Sigma _{\rm nl}} \langle \alpha \rangle |$|	max |Δc〈α〉|
MD-Patchy LRG (a)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0000	0.0006	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0010	0.0010	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0008	0.0007	–
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	0.9992	0.0004	0.0002
MD-Patchy LRG (b)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0020	0.0005	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0026	0.0010	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0017	0.0008	–
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0004	0.0007	0.0002
EZmock LRG (c)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0017	0.0005	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0019	0.0003	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0002	0.0007	–
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0018	0.0008	0.0004
EZmock ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0090	<10−4	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0087	0.0015	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0105	0.0016	–
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0103	<10−4	0.0003
EZmock LRG (c) × ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times \rbrace$|	1.0084	0.0004	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times \rbrace$|	1.0069	0.0001	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times ,\boldsymbol{\xi }_{\rm vv}^\times \rbrace$|	1.0067	0.0018	–
	|$\lbrace \boldsymbol{\xi }_{\rm comb}^\times \rbrace$|	1.0036	0.0006	0.0004

4.4 BAO fitting results

With the optimal void weights and our fiducial fitting ranges and priors detailed in Table 6, the fitting results from individual mock realizations are summarized in Table 9, for different data vectors. To examine the improvements on BAO constraints with voids, we measure the relative difference of σ_α as

$$\begin{eqnarray*} \delta \sigma _{\alpha , {\rm gg}} = \frac{ \sigma _{\alpha , {\rm gg}} - \sigma _\alpha }{\sigma _{\alpha , {\rm gg}} } , \end{eqnarray*}$$

(38)

where σ_{α, gg} indicates the fitted errors from galaxies alone. In addition, the best-fitting α parameters of each individual mock realization, as well as the fitted errors, are shown in Figs 12 and 13, respectively. In general, the dispersions of best-fitting α are consistent with the median fitted errors. It suggests that the statistical errors estimated from the fitted posteriors are unbiased.

Figure 13.

Distributions of the fitted error on α. Orange points indicate results from individual mocks, with different data vectors (see Fig. 11), and red histograms show their distributions. Purple plus denote measurements from fits to the mean 2PCFs of all mocks. Circles and dashed lines in blue indicate results from SDSS data.

Table 9.

A summary of BAO fitting results for the approximate mocks used in this work. m₁ and m₂ are the factors for correcting fitted parameter covariances, and the scatter of the best-fitting parameters of mocks, respectively (see Section 3.4.4). The symbol 〈 · 〉 denotes the median value of results from all 1000 individual mocks. Δα_{1σ −} and Δα_{1σ +} indicate the 16th and 84th percentile of the best-fitting α of individual mocks, and σ_α denotes the fitted error of α. σ_{α, syst} is the systematic error on α, obtained from the biases of α shown in Tables 7 and 8, which are added in quadrature. δσ_{α, gg} indicates the relative difference of the fitted error, compared to the result from galaxy–galaxy 2PCF of each mock [see equation (38)], and η_improve denotes the proportion of mocks with δσ_{α, gg} > 0, i.e. with smaller σ_α when including contributions from voids.

Sample	Data vector	m1	m2	|$\langle \alpha \rangle _{-\Delta \alpha _{1\sigma -}}^{+{\Delta \alpha _{1\sigma +}}}$|	〈σα〉	〈χ2〉/degree of freedom	σα, syst	〈δσα, gg〉	ηimprove
MD-Patchy LRG (a)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0173	1.0413	|$1.0000_{-0.0118}^{+0.0117}$|	0.0123	16.4/(22 − 5) = 0.96	0.0009	–	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0342	1.0807	|$1.0010_{-0.0110}^{+0.0100}$|	0.0109	32.1/(42 − 9) = 0.97	0.0014	11.7 %	97.6 %
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0494	1.1177	|$1.0008_{-0.0105}^{+0.0110}$|	0.0108	45.0/(60 − 13) = 0.96	0.0017	12.4 %	97.3 %
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0122	1.0318	|$0.9992_{-0.0106}^{+0.0107}$|	0.0107	11.4/(18 − 6) = 0.95	0.0005	13.7 %	98.2 %
MD-Patchy LRG (b)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0173	1.0413	|$1.0020_{-0.0107}^{+0.0114}$|	0.0109	16.3/(22 − 5) = 0.96	0.0021	–	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0342	1.0807	|$1.0026_{-0.0099}^{+0.0095}$|	0.0101	31.9/(42 − 9) = 0.97	0.0023	7.7 %	90.6 %
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0494	1.1177	|$1.0017_{-0.0098}^{+0.0102}$|	0.0100	45.0/(60 − 13) = 0.96	0.0020	8.0 %	91.8 %
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0122	1.0318	|$1.0004_{-0.0095}^{+0.0102}$|	0.0099	11.6/(18 − 6) = 0.97	0.0013	10.1 %	94.8 %
EZmock LRG (c)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0173	1.0413	|$1.0017_{-0.0145}^{+0.0159}$|	0.0149	16.4/(22 − 5) = 0.96	0.0013	–	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0342	1.0807	|$1.0019_{-0.0145}^{+0.0144}$|	0.0142	32.6/(42 − 9) = 0.99	0.0013	4.8 %	80.8 %
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0494	1.1177	|$1.0002_{-0.0107}^{+0.0112}$|	0.0107	46.3/(60 − 13) = 0.99	0.0024	27.7 %	100 %
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0122	1.0318	|$1.0018_{-0.0139}^{+0.0139}$|	0.0136	11.8/(18 − 6) = 0.98	0.0027	10.2 %	89.4 %
EZmock ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	1.0173	1.0413	|$1.0090_{-0.0364}^{+0.0440}$|	0.0409	16.7/(22 − 5) = 0.98	0.0057	–	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv}\rbrace$|	1.0342	1.0807	|$1.0087_{-0.0359}^{+0.0419}$|	0.0376	31.8/(42 − 9) = 0.96	0.0059	6.1 %	66.1 %
	|$\lbrace \boldsymbol{\xi }_{\rm gg},\boldsymbol{\xi }_{\rm gv},\boldsymbol{\xi }_{\rm vv}\rbrace$|	1.0494	1.1177	|$1.0105_{-0.0361}^{+0.0430}$|	0.0369	45.3/(60 − 13) = 0.96	0.0065	8.1 %	72.1 %
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	1.0122	1.0318	|$1.0103_{-0.0389}^{+0.0393}$|	0.0393	11.9/(18 − 6) = 0.99	0.0069	7.8 %	62.9 %
EZmock LRG (c) × ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times \rbrace$|	1.0173	1.0413	|$1.0084_{-0.0389}^{+0.0377}$|	0.0387	16.2/(22 − 5) = 0.95	0.0022	–	–
	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times \rbrace$|	1.0516	1.1224	|$1.0069_{-0.0359}^{+0.0376}$|	0.0360	46.3/(62 − 13) = 0.95	0.0024	6.8 %	74.3 %
	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times ,\boldsymbol{\xi }_{\rm gv}^\times ,\boldsymbol{\xi }_{\rm vg}^\times ,\boldsymbol{\xi }_{\rm vv}^\times \rbrace$|	1.0674	1.1616	|$1.0067_{-0.0342}^{+0.0378}$|	0.0351	59.0/(80 − 17) = 0.94	0.0046	8.7 %	75.8 %
	|$\lbrace \boldsymbol{\xi }_{\rm comb}^\times \rbrace$|	1.0122	1.0318	|$1.0036_{-0.0361}^{+0.0386}$|	0.0376	11.6/(18 − 6) = 0.97	0.0010	6.3 %	61.1 %

We sum the biases of α from the fits with different fitting ranges and parameter priors in quadrature to assess the systematic errors on the α measurements. Compared to the magnitude of statistical errors, systematic errors on the BAO measurements are typically fairly small. When including clustering measurements with voids, either by fitting simultaneously multiple 2PCFs, or by measuring BAO from the combined correlation function of galaxies and voids, the fitted error is smaller than that from galaxy–galaxy correlation alone, for the majority of the mock realizations. The typical improvement on the statistical uncertainty of BAO peak position is around 5–10 per cent. This confirms that as tracers of underdensities, the clustering of voids bring additional information compared to that of galaxies. Since our voids are defined by quadruples of galaxies, the information should be encoded in the galaxy distribution already, but exacted only from high-order statistics such as four-point correlations. Thus, the relative improvement from voids can be larger without BAO reconstruction, which transfers BAO information from higher-point statistics into two-point clustering (Schmittfull et al. 2015). Nevertheless, the joint BAO constraint after reconstruction should still be better, as long as the galaxy density field is not fully Gaussianized after reconstruction. This has been confirmed by the pre- and post-reconstruction measurements in Forero-Sánchez et al. (2021).

Furthermore, the improvement is relatively larger and detected from more mock realizations, when the BAO fitting uncertainty is smaller. This is consistent with the studies in Zhao et al. (2020), and suggests that the contribution of voids will be more promising for future surveys with larger volumes and smaller statistical errors. By comparing results from joint fits to the stacked data vectors with different clustering measurements, we conclude that the improvement of the BAO constraint is dominated by the contribution of galaxy–void cross-correlation. This is not surprising given the BAO signal-to-noise analysis in Section 3.2.3. However, for the ‘LRG (c)’ sample, an over 20 per cent improvement is observed when including voids. This is possibly due to the small fitted errors of void auto correlations, which are actually smaller than those of both galaxy–galaxy and galaxy–void correlations on average (see Table D1). We have examined different parameter priors for both the void 2PCFs of individual mocks, as well as the mean result from all mocks, but no obvious problem is found (see Appendix  D). In particular, for the mean void 2PCF of mocks, the fitted error with our fiducial priors is smaller than that with sufficiently large prior ranges, and agrees with the distribution of results from individual mocks. It implies that a prior is needed even for the fit to the mean void 2PCFs of the ‘LRG (c)’ mocks.

For the combined correlation function, however, the relative BAO precision improvement for ‘LRG (c)’ is around 10 per cent. This is more consistent with the results of the other samples. Actually, results from the combined 2PCFs are in good agreements with those from fits to multiple clustering measurements for all the other samples. This suggests that results from the combined 2PCFs are more robust. Moreover, these results benefit from the smaller lengths of data vector, compared to those of fits to the stacked 2PCFs. In this case, the correction factors for the uncertainties of the covariance matrices estimated using mocks are smaller, and statistical errors from the fits are more reliable. Therefore, we choose the multitracer fitting method based on the combined correlation functions with void weights, for BAO measurements and cosmological parameter constraints with the SDSS data, as is detailed in Section 5.

4.5 Cosmological parameter measurements

In order to perform cosmological analysis with the BAO fitting results, we convert the best-fitting α of individual mocks to |$D_{_{\rm V}} / r_{\rm d}$| measurements, and plot the correlation coefficients between different samples in Fig. 14. The results from galaxy–galaxy correlations and the combined 2PCFs are almost identical. This is consistent with the cross covariances of 2PCFs for different samples in Appendix  A, and can be explained by the fact that adding galaxy–void and void–void correlations does not bias the α measurement, but only reduces the fitted error. The correlation coefficient matrices can be divided into two blocks, one for the BOSS LRG samples that are below redshift 0.6, and the other for galaxies with z > 0.6. We do not consider the cross-correlations between these two populations, as they do not overlap in the 3D volume. Since the ‘LRG (a)’ and ‘LRG (b)’ samples share the same galaxies in the redshift range 0.4 < z < 0.5, their cross coefficient is relatively large. In contrast, the coefficient between the ‘LRG (c)’ and ‘ELG’ samples is small, as these two samples do not contain the same objects. Actually, these two samples are only correlated through their common volume, which is only a small proportion of the LRG volume, but involves the majority of the ELGs. This explains why the |$D_{_{\rm V}} / r_{\rm d}$| measured from the cross-correlations between ‘LRG (c)’ and ‘ELG’ is more correlated with the ‘ELG’ sample.

Figure 14.

Correlation coefficients of |$D_{_{\rm V}} / r_{\rm d}$| for different samples, based on the α measurements from 1000 individual mock realizations. ‘CROSS’ denotes results from the cross-correlations between ‘LRG (c)’ and ‘ELG’. ξ_gg and ξ_comb indicate results from galaxy–galaxy correlations and the combined 2PCFs, respectively.

We construct also the covariance matrices of |$D_{_{\rm V}}/r_{\rm d}$| from individual mocks, and rescaled them by the factor m₂ given by equation (33). The covariances for different samples are then used together with the median distance measurements of all mocks, for constraining H₀ and Ω_m in a flat-ΛCDM cosmological model. Apart from comparing measurements with galaxies alone and a combination of galaxies and voids, we examine also the galaxy-only results with and without the cross-correlations between the ‘LRG (c)’ and ‘ELG’ samples, to investigate improvements from multitracer analysis with different types of galaxies residing in the same cosmic volume. The constraints from combinations of BBN and BAO measurements are shown in Fig. 15 and Table 10. In all cases, the ‘true’ cosmological parameters of the mocks are well recovered. When including only the contribution from the cross-correlation between LRGs and ELGs, the uncertainties of both parameters are reduced by around 6 per cent, and a further 10 per cent improvement is achieved after taking into account the clustering of voids.

Figure 15.

Constraints on H₀ and Ω_m in a flat-ΛCDM model, obtained from the combination of BBN and different BAO measurements from mocks. |$\xi _{\rm gg}^{\rm auto}$| indicates results from galaxy auto correlation only, while |$\xi _{\rm gg}^{\rm all}$| includes also the contribution of galaxy cross-correlation between the ‘LRG (c)’ and ‘ELG’ samples. |$\xi _{\rm comb}^{\rm all}$| denotes BAO measurements from all samples, with both galaxies and voids. Dotted lines indicate the cosmological parameters for constructing the mocks. The contours show the 68 per cent and 95 per cent confidence intervals.

Table 10.

Best-fitting cosmological parameters and 1σ confidence intervals for the mocks, drawn from the posterior distributions in Fig. 15. The ‘true’ cosmological parameters denote the ones used for mock construction.

Parameter	‘True’	BAO (⁠|$\xi _{\rm gg}^{\rm auto}$|⁠)	BAO (⁠|$\xi _{\rm gg}^{\rm all}$|⁠)	BAO (ξcomb)
H0	67.77	|$68.7_{-5.1}^{+3.1}$|	|$68.7_{-4.7}^{+3.0}$|	|$68.0_{-4.1}^{+2.9}$|
Ωm	0.307115	|$0.330_{-0.095}^{+0.069}$|	|$0.329_{-0.087}^{+0.065}$|	|$0.313_{-0.076}^{+0.061}$|

The results confirm that the multitracer BAO analysis approach is promising for extracting more information from the cosmic density field, compared to the results with matter tracers alone. In addition, the unbiased cosmological parameter constraints with the mocks validate the robustness of our methodology, which can then be applied to the observational data.

5 RESULTS

As is discussed in the previous section, the multitracer BAO fitting approach based on the combined correlation functions is preferred over the joint fit to the stacked 2PCFs of multiple tracers. We then perform BAO fits with the combined 2PCFs of galaxies and voids from the SDSS data, with the fiducial BAO fitting scheme shown in bold in Table 6, and compare the results to those from galaxy correlations alone. In particular, when fitting to the galaxy–galaxy correlations, we use the Σ_nl values obtained from N-body simulations. The fitting results are converted to distance measurements, which are finally used for cosmological parameter constraints.

5.1 Distance measurements

With the fit settings listed in Table 6, we run BAO fits for both the galaxy–galaxy correlation (ξ_gg) and combined correlation functions with galaxies and voids (ξ_comb) for all the SDSS samples studied in this work, with covariance matrices estimated using the corresponding approximate mocks. The best-fitting BAO model curves agree well with the observational data in the fitting ranges, as shown in Fig. 16, with the corresponding posterior distributions of the α parameters illustrated in Fig. 17. The best-fitting α parameters, their statistical uncertainties drawn from the posterior distributions, as well as the reduced chi-squared values, are listed in Table 11. For all the samples, the α measurements from ξ_gg and ξ_comb agree within 0.5σ statistical uncertainties. Besides, the fitted errors for the combined correlation functions are always smaller than those from galaxies alone, with the relative differences ranging from 5 per cent to 15 per cent.

Figure 16.

Galaxy–galaxy correlations and combined 2PCFs measured from SDSS data, as well as the best-fitting model curves with and without BAO peaks. In particular, the combined correlations are computed with the optimal void weights measured from mocks. Shadowed areas indicate regions outside the fitting ranges.

Figure 17.

Marginalized posterior distributions of α for different SDSS samples. Orange lines show results from galaxies only, while purple lines indicate measurements from the combined 2PCFs of galaxies and voids.

Table 11.

BAO fitting results for the SDSS data. The best-fitting value and lower and upper 1σ confidence limits are evaluated from the 50th, 16th, and 84th percentiles of the marginalized posterior distribution of α. The systematic error σ_{α, syst} is estimated using mocks, as is shown in Table 9. (χ²/degree of freedom) indicates the reduced chi-squared obtained from the model with BAO peaks, and |$\Delta \chi ^2 = \chi _{\rm no\, BAO}^2 - \chi _{\rm BAO}^2$|⁠, where |$\chi _{\rm no\, BAO}^2$| and |$\chi _{\rm BAO}^2$| denote the best-fitting chi-squared values from the model without and with BAO peaks, respectively.

Sample	Data	αfit	σα, syst	|$\dfrac{\chi ^2 }{\rm d.o.f.}$|	Δχ2
	vector
LRG (a)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	|$0.9992_{-0.0097}^{+0.0100}$|	0.0009	1.51	55.9
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	|$0.9952_{-0.0086}^{+0.0083}$|	0.0005	1.59	62.8
LRG (b)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	|$0.9928_{-0.0089}^{+0.0088}$|	0.0021	1.18	57.2
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	|$0.9927_{-0.0076}^{+0.0076}$|	0.0013	0.46	70.7
LRG (c)	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	|$0.9982_{-0.0119}^{+0.0116}$|	0.0013	1.46	37.9
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	|$0.9934_{-0.0109}^{+0.0110}$|	0.0027	1.26	41.0
ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}\rbrace$|	|$0.9825_{-0.0318}^{+0.0321}$|	0.0057	1.34	2.2
	|$\lbrace \boldsymbol{\xi }_{\rm comb}\rbrace$|	|$0.9867_{-0.0236}^{+0.0310}$|	0.0069	0.98	5.7
LRG (c) × ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times \rbrace$|	|$0.9565_{-0.0246}^{+0.0254}$|	0.0022	1.10	14.3
	|$\lbrace \boldsymbol{\xi }_{\rm comb}^\times \rbrace$|	|$0.9539_{-0.0212}^{+0.0228}$|	0.0010	0.56	11.6

The improvements on the BAO measurements are significant, compared to the systematic errors estimated using approximate mocks, as well as uncertainties on the fitted errors, which are generally less than 5 per cent (see Appendix  D). Fig. 16 shows also the best-fitting models without BAO peaks, obtained by replacing P_lin(k) with P_{lin, nw}(k) in equation (21). The smoothness of best-fitting ‘non-wiggle’ curves confirm that our BAO models for both the galaxy–galaxy and combined 2PCFs are able to describe well the BAO peak, without introducing additional patterns on the same scale. The differences of chi-squared values, for models without and with BAO peaks, are presented in Table 11. These results suggest that the BAO significances of the combined correlation functions are generally higher than the ones of ξ_gg, despite the smaller number of ξ_comb bins.

The fitting results of the SDSS data are also shown in Figs 12 and 13, in comparison with results from the approximate mocks. The best-fitting α values from the data and mocks are generally in good agreements. However, the fitted errors of the data, computed as the mean of the lower and upper 1σ confidence limits, are always smaller than the median BAO errors of mocks. This is consistent with previous studies (e.g. Vargas-Magaña et al. 2018; Zhao et al. 2020; Bautista et al. 2021), and can be explained by the fact that the approximate mocks generally overestimate the nonlinear damping of the BAO peak, which magnifies the uncertainty of BAO position. We have further checked results for different Galactic caps, and good consistencies are found in general (see Appendix  E).

Given the |$D_{_{\rm V}} / r_{\rm d}$| value in our fiducial cosmology at different redshifts, we convert the fitting results of α from the SDSS data to cosmological distance measurements, and present the results in Table 12. Since the upper and lower 1σ errors of α are similar for almost all cases, we quote only symmetric errors of |$D_{_{\rm V}} / r_{\rm d}$|⁠. Moreover, the statistical errors are rescaled by |$\sqrt{m_1}$|⁠, following equation (32), and then combined with systematic errors by summing in quadrature. The cross covariances between different samples are estimated using the correlation coefficients measured from 1000 mock realizations (see Fig. 14), and are presented in Appendix  F.

Table 12.

Distance measurements from the BAO fitting results of the SDSS data at different effective redshifts. Systematic errors estimated using mocks are also included, for both results from galaxies alone (ξ_gg), and measurements from a joint analysis with galaxies and voids (ξ_comb).

Sample	zeff	ξgg	ξcomb
		|$D_{_{\rm V}} / r_{\rm d}$|	|$D_{_{\rm V}} / r_{\rm d}$|	Precision
LRG (a)	0.38	9.98 ± 0.10	9.94 ± 0.08	0.9 %
LRG (b)	0.51	12.67 ± 0.12	12.67 ± 0.10	0.8 %
LRG (c)	0.70	16.25 ± 0.19	16.17 ± 0.18	1.1 %
LRG (c) × ELG	0.77	16.67 ± 0.44	16.63 ± 0.39	2.3 %
ELG	0.85	18.33 ± 0.61	18.41 ± 0.53	2.9 %

We then compare our distance measurements from the combined correlation functions with those from previous studies – that are based on data from different spectroscopic surveys – in Fig. 18. In particular, the measurements are collected from BAO-only analysis of 6dFGS (Beutler et al. 2011), WiggleZ (Kazin et al. 2014), SDSS MGS (Ross et al. 2015), BOSS LRGs (Alam et al. 2017), eBOSS LRGs (Gil-Marín et al. 2020; Bautista et al. 2021), eBOSS ELGs in configuration space (Raichoor et al. 2021), and Fourier space (de Mattia et al. 2021), eBOSS QSOs (Neveux et al. 2020; Hou et al. 2021), as well as eBOSS Lyα forests (du Mas des Bourboux et al. 2020). Here, the consensus results based on configuration and Fourier space measurements are used whenever possible, and the Lyα forest measurement is based on the combination of Lyα auto correlations and Lyα–QSO cross-correlations. Previous BOSS/eBOSS results that are given as |$D_{_{\rm M}}/r_{\rm d}$| and |$D_{_{\rm H}}/r_{\rm d}$| measurements, are combined into |$D_{_{\rm V}}/r_{\rm d}$| following equation (19), with their cross covariances taken into account. The distance measurements are then normalized by predictions of the best-fitting base-ΛCDM parameters from Planck Collaboration (2020), for better comparisons of the sizes of error bars. Besides, measurements at the same effective redshifts are added small vertical offsets for illustration purposes.

Figure 18.

BAO distance measurements from different spectroscopic data, normalized by predictions of the Planck Collaboration (2020) best-fitting base-ΛCDM parameters (Ω_m = 0.3158, h = 0.6732). The results are all from BAO-only analysis.

In general, the results in this work agree well with previous BOSS/eBOSS analyses, apart from the smaller uncertainties due to the contributions of voids. One can also see that the central values of our measurements are always smaller than the predictions of Planck results, though the differences are not significant. In particular, the result from the cross-correlations between ‘LRG (c)’ and ‘ELG’ shows the largest tension, which is around 2|$\, \sigma$|⁠. We then study the implications on cosmological parameters in the following sections.

5.2 Hubble parameter and density parameters

To highlight the contribution of the multitracer approach, we compare firstly the cosmological parameters constrained using the distances measured from the LRG and ELG samples in this work, to the corresponding measurements used in eBOSS21. In addition, we combine the |$D_{_{\rm M}} / r_{\rm d}$| and |$D_{_{\rm H}} / r_{\rm d}$| measured from BOSS/eBOSS LRGs into |$D_{_{\rm V}} / r_{\rm d}$|⁠, to have 1D BAO results for fair comparisons. The posterior distributions of H₀, Ω_m, and |$\Omega _\Lambda h^2$| from joint BBN and BAO analyses in a flat-ΛCDM cosmology are shown in Fig. 19. Note that throughout this work, the contours show the 68 per cent and 95 per cent confidence intervals. It can be seen that for the 1D BAO constraints based on |$D_{_{\rm V}} / r_{\rm d}$|⁠, the multitracer results with cosmic voids improve significantly the precision of all parameters. However, the uncertainties of parameters measured from multitracer 1D BAOs can be larger than the ones from 2D BAOs of galaxies alone. For instance, the marginalized 1σ error of the Ω_m parameter from the multitracer 1D BAOs is around 18 per cent smaller than the measurements from the 1D BAOs of galaxies alone, but 19 per cent larger than the corresponding galaxy 2D BAO results. This suggests that the decomposed |$D_{_{\rm M}} / r_{\rm d}$| and |$D_{_{\rm H}} / r_{\rm d}$| parameters encode more cosmological information than the spherical averaged |$D_{_{\rm V}} / r_{\rm d}$|⁠.

Figure 19.

Constraints on H₀, Ω_m, and |$\Omega _\Lambda h^2$| in a flat-ΛCDM model, obtained from the combination of BBN and different BAO measurements. Orange contours and lines show the results with |$D_{_{\rm M}} / r_{\rm d}$| and |$D_{_{\rm H}} / r_{\rm d}$| measured from LRGs and ELGs in eBOSS21. Red-dashed lines indicate constraints with the corresponding |$D_{_{\rm V}} / r_{\rm d}$|⁠, obtained by combining the 2D distance measurements, with their cross covariances taken into account. Black solid lines indicate results with |$D_{_{\rm V}} / r_{\rm d}$| measured in a multitracer manner, using the combined correlation functions of galaxies and voids, including the cross-correlations between LRGs and ELGs.

Nevertheless, it is still possible that the multitracer 1D BAOs yields better constraints for certain parameters, compared to the results from 2D BAOs of galaxies alone. As an example, the anisotropic BAO do not constrain Ω_Λh² better than the isotropic BAO in the flat-ΛCDM framework at low redshifts, since both |$D_{_{\rm M}} / r_{\rm d}$| and |$D_{_{\rm H}} / r_{\rm d}$| depend on this parameter in approximately the same way, in the dark-energy-dominated era. Actually, we observe similar |$\Omega _\Lambda h^2$| uncertainties for both 1D and 2D galaxy BAOs in Fig. 19. In contrast, a 16 per cent reduction of the 1σ error is observed for the 1D BAO with voids. This implies that by combining the multitracer 1D BAOs with different data sets, the cosmological constraints can be tighter than those from the galaxy 2D BAOs, due to different parameter degeneracies. We start with combining the multitracer BAOs measured in this work, with the BAO-only results from the rest of the SDSS large-scale structure samples, including MGS (Ross et al. 2015), QSOs (Neveux et al. 2020; Hou et al. 2021), and Lyα forests (du Mas des Bourboux et al. 2020). The resulting BAO measurements are dubbed ‘mtBAO’ henceforth. It is worth noting that the same data sets are used for cosmological parameter constraints with BAO in eBOSS21, which serves as our baseline for examining the contributions from cosmic voids.

The joint BBN and BAO constraints in a flat-ΛCDM cosmology are shown in Fig. 20, for both ‘mtBAO’ and eBOSS21 BAO in different redshift ranges. In general, the results from ‘mtBAO’ and eBOSS21 BAO are consistent. The left-hand panel shows that the |$\Omega _\Lambda h^2$| constraint is dominated by the low-redshift (z < 1) BAO measurements, while BAOs of tracers with effective redshifts above one contribute mainly to the Ω_mh² measurement. As the result, with ‘mtBAO’ the 1σ error of |$\Omega _\Lambda h^2$| is reduced by ∼17 per cent, compared to the results with eBOSS21 BAO, while the uncertainties of the Ω_mh² parameter are similar for the two cases. The same posterior distributions projected on the H₀–Ω_m plane are shown in the right panel of Fig. 20, with the marginalized results listed in Table 13. Though the 1D ‘mtBAO’ at z < 1 yields worse constraints for both parameters than the ones from the 2D eBOSS21 BAO, the overall 1σ uncertainties of H₀ and Ω_m, with the QSOs and Lyα forests taken into account, are both reduced by around 6 per cent.

Figure 20.

Constraints on Ω_Λh², Ω_mh², H₀, and Ω_m in a flat-ΛCDM model, obtained from the combination of BBN and different BAO measurements. Orange, red, and blue contours show results taken from eBOSS21, for SDSS 2D BAO measurements in different redshift bins. Dashed and solid lines indicate the multitracer BAO measurements of galaxies and voids in this work, together with measurements from the other SDSS tracers, including MGS, QSOs, and Lyα forests. Grey bands refer to the (⁠|$68\, {{\ \rm per\ cent}}, 95\, {{\ \rm per\ cent}}$|⁠) H₀ measurements from the distance ladder method (Riess et al. 2019).

Table 13.

Marginalized best-fitting values and 68 per cent confidence limits of cosmological parameters in the flat-ΛCDM model and its two one-parameter extensions, obtained using different probes, including Planck CMB temperature and polarization data (Planck Collaboration 2020), Pantheon SNe Ia sample (Scolnic et al. 2018), and distance ladder measurements (Riess et al. 2019). ‘BAO’ denotes the SDSS BAO results used in eBOSS21, where the BOSS/eBOSS LRG BAOs are measured in 2D. In contrast, ‘mtBAO‘ indicates BAO measurements from the same SDSS data, but with the BOSS/eBOSS LRG and ELG results replaced by the multitracer 1D BAO measurements in this work.

Model	Probe	H0 (⁠|${\rm km}\, {\rm s}^{-1}\, {\rm Mpc}^{-1}$|⁠)	Ωm	|$\Omega _\Lambda h^2$|	Ωk	w
ΛCDM	BBN + BAO	67.35 ± 0.98	0.299 ± 0.016	0.3179 ± 0.0094	–	–
	BBN + mtBAO	67.58 ± 0.91	0.290 ± 0.015	0.3241 ± 0.0079	–	–
	CMB	67.29 ± 0.61	0.3164 ± 0.0084	0.3096 ± 0.0094	–	–
	CMB + BAO	67.60 ± 0.43	0.3119 ± 0.0058	0.3145 ± 0.0066	–	–
	CMB + mtBAO	67.96 ± 0.39	0.3070 ± 0.0051	0.3201 ± 0.0060	–	–
	Distance ladder	74.0 ± 1.4	–	–	–	–
	Distance Ladder + BAO	74.0 ± 1.4	0.299 ± 0.016	0.384 ± 0.017	–	–
	Distance Ladder + mtBAO	74.0 ± 1.4	0.290 ± 0.015	0.390 ± 0.017	–	–
oΛCDM	BAO	–	0.285 ± 0.023	|$0.35_{-0.20}^{+0.11}$|	|$0.079_{-0.10}^{+0.083}$|	–
	mtBAO	–	0.284 ± 0.022	|$0.38_{-0.21}^{+0.12}$|	|$0.04_{-0.12}^{+0.095}$|	–
	CMB	54.5 ± 3.6	|$0.483_{-0.069}^{+0.055}$|	|$0.169_{-0.041}^{+0.030}$|	|$-0.044_{-0.014}^{+0.019}$|	–
	SN	–	0.317 ± 0.071	–	−0.05 ± 0.18	–
	SN + BAO	–	0.289 ± 0.021	|$0.38_{-0.21}^{+0.11}$|	0.037 ± 0.069	–
	SN + mtBAO	–	0.288 ± 0.020	|$0.39_{-0.21}^{+0.12}$|	0.020 ± 0.070	–
wCDM	BAO	–	|$0.271_{-0.017}^{+0.038}$|	|$0.38_{-0.24}^{+0.11}$|	–	−0.69 ± 0.15
	mtBAO	–	|$0.256_{-0.016}^{+0.048}$|	|$0.39_{-0.25}^{+0.12}$|	–	−0.65 ± 0.16
	CMB	–	|$0.199_{-0.057}^{+0.022}$|	|$0.61_{-0.11}^{+0.23}$|	–	|$-1.58_{-0.35}^{+0.16}$|
	SN	–	|$0.309_{-0.056}^{+0.088}$|	–	–	|$-1.07_{-0.18}^{+0.24}$|
	SN + BAO	–	0.292 ± 0.015	|$0.39_{-0.23}^{+0.11}$|	–	−0.950 ± 0.061
	SN + mtBAO	–	0.286 ± 0.014	|$0.40_{-0.23}^{+0.12}$|	–	−0.943 ± 0.060

The H₀ value measured from BBN and BAO is consistent with the prediction of Planck CMB temperature and polarization data (Planck Collaboration 2020), but significantly smaller than the local distance ladder measurement (Riess et al. 2019). This is known as the ‘H₀ tension’. To understand better the contributions of multitracer BAOs to the H₀ measurements, we further examine the constraints on H₀ and the sound horizon r_d from different probes, as well as their combinations. The 2D confidence intervals and marginalized statistics are shown in Fig. 21 and Table 13, respectively. It can be seen that BAO only constrains |$(h\, r_{\rm d})$|⁠, or r_d expressed in the unit of |$h^{-1}\, {\rm Mpc}$|⁠. Thus, it does not improve the precision of H₀ measurement when combined with distance ladder results. Nevertheless, the constraint on |$(h\, r_{\rm d})$| can be translated to that of H₀, when r_d measurements are available. This can be achieved with the CMB data, or through the ω_b value constrained by BBN. Actually, the ω_b measured from BBN is consistent with that from the Planck CMB data (Cooke et al. 2018). Therefore, they produce consistent H₀ constraints when combined with BAO. Compared to the results with eBOSS21 BAO, the 1σ uncertainty of |$(h\, r_{\rm d})$| measured from ‘mtBAO’ is around 12 per cent smaller. This corresponds to 6 per cent or 10 per cent reduction of the H₀ uncertainty, when combined with BBN or CMB data, respectively.

Figure 21.

Constraints on |$(h\, r_{\rm d})$| and H₀ in a flat-ΛCDM model, for SDSS BAO measurements in eBOSS21 (‘BAO’), and the ones with multitracer 1D BAO from LRGs and ELGs measured in this work (‘mtBAO’, see also Fig. 20). The results are also combined with BBN, the Planck CMB temperature and polarization data, or the distance ladder measurements.

5.3 Curvature and dark energy

We further examine the role of multitracer BAO measurements on two single-parameter extensions to the flat-ΛCDM model, i.e. oΛCDM and wCDM, which are the cosmologies with free curvature and dark energy equation of state, respectively (see Section 3.5). Following eBOSS21, we include also the Planck CMB temperature and polarization data (Planck Collaboration 2020) and the Pantheon SNe Ia sample (Scolnic et al. 2018). The results from single probes, as well as the combination of BAO and SNe Ia data, are shown in Figs 22 and 23, with the marginalized statistics listed in Table 13. Among these measurements, BAO provides the tightest single-probe constraints on the extension parameters Ω_k and w, which are also less coupled with Ω_m compared to the results from CMB and SNe Ia data. For both cosmological models, the posterior contours from BAOs are in around 2σ tensions with those from CMB alone. Therefore, we combine only the BAO and SNe Ia measurements, which are in better agreements. With the combined constraints, a flat-ΛCDM model is favoured for both extended models.

Figure 22.

Constraints on |$\Omega _\Lambda$| and Ω_m in a ΛCDM model with free curvature (oΛCDM), for Planck CMB data, Pantheon SNe Ia sample, and the same SDSS BAO measurements as in Fig. 21. The dash–dotted line indicates a flat-ΛCDM model.

Figure 23.

Same as Fig. 22, but for w and Ω_m constraints in a flat-wCDM cosmology model.

For the BAO-only constraints on the Ω_k and w parameters, the 1σ errors from the multitracer 1D BAOs are both over 10 per cent larger than those from the eBOSS21 2D BAOs. It indicates that the anisotropic BAO measurements are crucial for the extended cosmological parameters. Nevertheless, after combining with the SNe Ia data, the precisions of these parameters measured from ‘mtBAO’ and eBOSS21 BAO are comparable. In terms of the Ω_m parameter, for both models the 1σ uncertainties from ‘mtBAO’ are ∼6 per cent smaller the ones from eBOSS21 BAO, when the SNe Ia data are involved.

6 CONCLUSIONS

We have performed a multitracer BAO analysis with different types of galaxies, as well as the corresponding cosmic voids identified by the dive algorithm (Zhao et al. 2016), based on the SDSS BOSS DR12 and eBOSS DR16 data. The data sets include over 1.6 million LRGs spanning in the redshift range of 0.2 < z < 1.0 (Reid et al. 2016; Ross et al. 2020), and 0.17 million star-forming ELGs at 0.6 < z < 1.1 (Raichoor et al. 2021). In addition, one thousand realizations of multitracer approximate mock catalogues are used for each galaxy sample (Kitaura et al. 2016b; Zhao et al. 2021a), for method validations and covariance matrix estimations. Since the footprints of the two types of galaxies overlap, we take into account also their cross-correlation functions.

The void catalogues are constructed upon galaxy samples with the BAO reconstruction technique applied, which improves the BAO peak significance by suppressing nonlinear damping effects through the reverse of bulk flow motions. To avoid contamination from voids-in-clouds which reside in overdensities (Liang et al. 2016; Zhao et al. 2016), we select voids by the radii and use only large spheres for clustering analysis. The radius thresholds are chosen to maximize the BAO SNR measured from mocks. As the result, the numbers of voids used in our studies are about twice the numbers of the corresponding galaxies.

To measure the BAO peak position with the contributions of both galaxies and voids, we have explored two different BAO fitting methods, including joint fits to galaxy–galaxy, galaxy–void, and void–void correlations with the same BAO dilation parameter α, as well as combining individual correlation functions for the fit, with a negative weight applied to voids (see Zhao et al. 2020). The optimal void weight is chosen to minimize the BAO fitting error for the mean combined 2PCFs of all mocks. We find no significant differences on the BAO measurements from both methods for most cases, and the systematic errors estimated using mocks are generally small compared to the corresponding statistical uncertainties. We then use the BAO fitting method with the combined 2PCFs for the SDSS data, due to the smaller length of the data vector, with which the covariance estimations are more accurate with a fixed number of mock realizations.

With the combined correlation functions of galaxies and voids, we observe on average around 10 per cent improvements on the precision of BAO measurements compared to the results from galaxies alone, for the mocks of all data sets. Besides, for the majority of the individual mock realizations, improvements from the clustering of voids are detected. It confirms that the clustering of underdensities encode additional information compared to that of galaxies alone. Since the DT voids are defined by galaxies in a model-independent manner, the extra information is possibly exacted from higher-point clustering of the galaxy sample.

Interestingly, there are more mocks with improvements when the BAO uncertainty of the galaxy sample is smaller. This is consistent with findings in Zhao et al. (2020), and suggests that the contribution from voids may be more promising for future surveys with larger data sets, such as DESI. In addition, cosmological fits in a flat-ΛCDM model recover well the H₀ and Ω_m parameters of the mocks. Both the cross-correlation between different types of galaxies, and the clustering of voids, are found to improve the precision of the constraints. This confirms that a multitracer BAO analysis reduces cosmic variances.

For the SDSS data, the statistical uncertainties of BAO measurements are reduced by around per cent 5 to 15 per cent for all the samples, when including the clustering statistics with voids. The corresponding distance measurements are in good agreements with the BAO analyses of previous studies using the same galaxy samples (e.g. Alam et al. 2017; Bautista et al. 2021; Raichoor et al. 2021). Although only isotropic BAO measurements are considered in this work, the cosmological constraints can be better than those from the anisotropic BAO analysis with galaxies alone for certain parameters, especially when the results are combined with different data sets, such as the Planck CMB data and the Pantheon SNe Ia sample. For instance, with the |$D_{_{\rm V}} / r_{\rm d}$| values measured in this work, the uncertainties of the H₀, Ω_m, and |$\Omega _\Lambda h^2$| parameters from a joint BBN and BAO analysis in the flat-ΛCDM cosmology, are reduced by around 6, 6, and 17 per cent respectively, compared to the results in eBOSS21. We foresee tighter cosmological constraints in the multi-tracer approach with anisotropic BAO measurements from cosmic voids, which may help in distinguishing different cosmological models. We leave relevant studies to a future work.

ACKNOWLEDGEMENTS

CZ, AV, DFS, AT, and JY acknowledge support from the Swiss National Science Foundation (SNF) grant 200020_175751. FSK acknowledges the grants SEV-2015-0548, RYC2015-18693, and AYA2017-89891-P. CT is supported by Tsinghua University and sino french CNRS-CAS international laboratories LIA Origins and FCPPL. GR acknowledges support from the National Research Foundation of Korea (NRF) through grant nos. 2017R1E1A1A01077508 and 2020R1A2C1005655 funded by the Korean Ministry of Education, Science, and Technology (MoEST).

Clustering measurements and BAO fits for this work were performed at the National Energy Research Scientific Computing Center (NERSC),²² U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

DATA AVAILABILITY

The SDSS galaxy data and the corresponding approximate mock catalogues used in this work can be downloaded from the SDSS Science Archive Server (see Section 2 for details). All software packages for the analysis, including the void finder, BAO fitter, and cosmological parameter sampler are all publicly available (see Section 3). In addition, the void catalogues can be obtained directly upon request to CZ.

Footnotes

1

https://github.com/cheng-zhao/DIVE

2

The BAO and RSD measurements are available at https://sdss.org/science/final-bao-and-rsd-measurements, and the cosmological results can be found in https://sdss.org/science/cosmology-results-from-eboss and https://svn.sdss.org/public/data/eboss/DR16cosmo/tags/v1_0_1/.

3

https://data.sdss.org/sas/dr12/boss/lss/

4

https://data.sdss.org/sas/dr16/eboss/lss/catalogs/DR16/

5

https://data.sdss.org/sas/dr12/boss/lss/dr12_multidark_patchy_mocks/

6

https://data.sdss.org/sas/dr17/eboss/lss/EZmocks/

7

https://github.com/mockFactory/make_survey

8

https://github.com/cheng-zhao/brickmask

9

https://github.com/julianbautista/eboss_clustering

10

https://github.com/seshnadathur/Revolver

11

https://www.cgal.org

12

https://github.com/cheng-zhao/FCFC

13

https://camb.info/

14

https://github.com/cheng-zhao/CosmoGAME

15

https://github.com/cheng-zhao/powspec

16

https://github.com/cheng-zhao/BAOflit

17

https://github.com/farhanferoz/MultiNest

18

https://github.com/JohannesBuchner/PyMultiNest

19

https://github.com/evamariam/CosmoMC_SDSS2020

20

http://cosmologist.info/cosmomc

21

https://github.com/cmbant/getdist

22

https://ror.org/05v3mvq14

APPENDIX A: COVARIANCES OF CLUSTERING MEASUREMENTS

Since voids are indirect tracers resolved from galaxy samples, the distributions of voids and galaxies are highly correlated. To quantify the correlations of their clustering measurements, we measure the correlation matrices of galaxy–galaxy (ξ_gg), galaxy–void (ξ_gv), and void-void (ξ_vv) correlations from the 1000 mock realizations for each SDSS sample, and plot the results in Fig. A1. Moreover, since the pair counts for these 2PCFs are used for the evaluations of combined correlation functions of galaxies and voids (ξ_comb; see Section 3.3), the correlation coefficients involving ξ_comb, which are computed with the optimal void weights found in Section 4.1, are also shown. In particular, we present results only in the separation range of |$(50,150)\, h^{-1}\, {\rm Mpc}$|⁠, which is close to the fiducial BAO fitting range used in this work (see Section 4.2).

Figure A1.

Correlation coefficients of galaxy–galaxy, galaxy–void, void–void, and the combined correlation functions in |$(50, 150)\, h^{-1}\, {\rm Mpc}$|⁠, measured from all mock realizations for each sample. Here, the combined correlation functions are computed with the optimal void weights found in Section 4.1.

As expected, the cross-correlation coefficients between different clustering measurements are generally significant. A typical diagonal correlation coefficient between ξ_gg and ξ_gv is −0.4. The negative value is due to the negative void biases (Zhao et al. 2016; Chuang et al. 2017). For the cross coefficients between ξ_gg and ξ_vv, the diagonal values are usually around 0.2. In addition, ξ_comb is strongly correlated with the individual 2PCFs, as its diagonal coefficients of the cross covariances with ξ_gg, ξ_gv, and ξ_vv are roughly 0.9, 0.7, and 0.4 respectively. These figures suggest that ξ_gg dominates ξ_comb, which is consistent with the small amplitude of the optimal void weights.

We further examine the cross covariances between clustering measurements of different samples. In particular, we consider the cross-correlations between ‘LRG (a)’ and ‘LRG (b)’, as they share the same galaxies and voids due to the overlap of their redshift ranges. Moreover, the cross covariances for ‘LRG (c)’, ‘ELG’, and their cross-correlations are taken into account as well. The results for ξ_gg and ξ_comb are shown in Fig. A2, where ξ_comb is again computed with the optimal weights found in Section 4.1. It can be seen that the correlation coefficients between ‘LRG (a)’ and ‘LRG (b)’ are relatively large. Besides, clustering statistics of the ‘ELG’ sample are strongly correlated with the cross-correlations between ‘LRG (c)’ and ‘ELG’. This can be explained by the fact that a large fraction of the ‘ELG’ sample is involved for the cross-correlations, while many of the LRGs are too far away from the ELG footprints for correlations on BAO scale (see Fig. 1).

Figure A2.

Correlation coefficients of correlation function measurements in |$(50, 150)\, h^{-1}\, {\rm Mpc}$|⁠, from different mock samples, for both galaxy 2PCFs (ξ_gg) and combined correlations with the optimal void weights (ξ_comb). ‘CROSS’ denotes the cross-correlations between ‘LRG (c)’ and ‘ELG’.

APPENDIX B: LEAST SQUARES FITTING AND CHI-SQUARED EVALUATION BASED ON QR DECOMPOSITION

Given a data vector |$\boldsymbol{y}$|⁠, the best-fitting parameters of the the model

$$\begin{eqnarray*} f(\boldsymbol{x}) = [ {\bf X} (\boldsymbol{x}) ]^{\sf T} \boldsymbol{a} \end{eqnarray*}$$

(B1)

can be solved linearly using the least squares method. Here, |$\boldsymbol{a}$| indicates the free parameters to be fitted, and |${\bf X} (\boldsymbol{x})$| denotes a set of arbitrary fixed functions of |$\boldsymbol{x}$|⁠. This is because the parameters that minimize the chi-squared function for model f,

$$\begin{eqnarray*} \chi _f^2 = (\boldsymbol{y} - {\bf X} \boldsymbol{a})^{\sf T} {\bf C}^{-1} (\boldsymbol{y} - {\bf X} \boldsymbol{a}) , \end{eqnarray*}$$

(B2)

are the solution of the following linear equation (e.g. Press et al. 2007)

$$\begin{eqnarray*} ({\bf X}^{\sf T} {\bf C}^{-1} {\bf X}) \cdot \boldsymbol{a} = {\bf X}^{\sf T} {\bf C}^{-1} \boldsymbol{y} , \end{eqnarray*}$$

(B3)

where C denotes the covariance matrix. Direct solution of this equation, as well as the inverse of the covariance matrix, can be numerically unstable. Moreover, since the chi-squared function (equation (25)) is typically evaluated again and again, with different parameters during a fitting process, solving equation (B3) directly can be computational expensive in practice, especially when the dimension of |$\boldsymbol{a}$| is high.

If the covariance matrix is obtained in the form of |${\bf C} = {\bf M}^{\sf T} {\bf M}$|⁠, however, the problem can be largely simplified. This is actually a common occurrence, as covariance matrices estimated using mocks and jackknife resampling can both be expressed in this way. The constant factors, such as the Hartlap correction (equation 30), can be seen as being absorbed by the M term here. This covariance generation matrix M, can be QR-decomposed as

$$\begin{eqnarray*} {\bf M} = {\bf Q R} , \end{eqnarray*}$$

(B4)

where Q is an orthogonal matrix (⁠|${\bf Q}^{\sf T} {\bf Q} = {\bf I}$|⁠, with I being the identity matrix), and R is an upper triangular matrix. Thus,

$$\begin{eqnarray*} {\bf C}^{-1} = ({\bf R}^{\sf T} {\bf Q}^{\sf T} {\bf Q R})^{-1} = {\bf R}^{-1} {\bf R}^{- {\sf T}} . \end{eqnarray*}$$

(B5)

equation (B3) can then be written as

$$\begin{eqnarray*} ({\bf R}^{- {\sf T}} {\bf X})^{\sf T} ({\bf R}^{- {\sf T}} {\bf X}) \cdot \boldsymbol{a} = ({\bf R}^{- {\sf T}} {\bf X})^{\sf T} {\bf R}^{- {\sf T}} \boldsymbol{y}. \end{eqnarray*}$$

(B6)

Let |${\bf Y} = {\bf R}^{- {\sf T}} {\bf X}$|⁠, it can be further QR-decomposed as Y = Q′U, where U is another upper triangular matrix. Then, equation (B6) can be simplified as

$$\begin{eqnarray*} ({\bf U}^{\sf T} {\bf U}) \cdot \boldsymbol{a} = {\bf Y}^{\sf T} {\bf R}^{- {\sf T}} \boldsymbol{y}. \end{eqnarray*}$$

(B7)

Let |${\bf Z} = {\bf U}^{-{\sf T}} {\bf Y}^{\sf T} {\bf R}^{- {\sf T}}$|⁠, we have finally

$$\begin{eqnarray*} {\bf U} \cdot \boldsymbol{a} = {\bf Z} \boldsymbol{y} . \end{eqnarray*}$$

(B8)

In practice, since M and X are known, the matrices U and Z can be pre-computed. Typically, |$\boldsymbol{y}$| is the fitting residual with the non-nuisance parameters. In this case, during the fitting procedure, one merely needs a matrix-vector production |$\boldsymbol{b} = {\bf Z} \boldsymbol{y}$|⁠, and solving the equation

$$\begin{eqnarray*} {\bf U} \cdot \boldsymbol{a} = \boldsymbol{b} . \end{eqnarray*}$$

(B9)

Since U is an upper triangular matrix, this equation can be solved efficiently with a backward substitution. It is worth noting that Y and Z can be evaluated using forward and backward substitutions as well, in a column-by-column manner, to avoid the inverse of upper and lower triangular matrices. Therefore, the whole least squares fitting process involves only QR decomposition, matrix production, as well as forward and backward substitutions. In particular, the two QR decompositions can both be computed in-place, as we do not use the orthogonal matrices explicitly. Thus, the algorithm here is highly numerically stable and efficient.

Furthermore, given the inverse covariance matrix expressed by equation (B5), the chi-squared function (equation 25) becomes

$$\begin{eqnarray*} \chi ^2 = \boldsymbol{d}^{\sf T} {\bf R}^{-1} {\bf R}^{- {\sf T}} \boldsymbol{d} = \Vert {\bf R}^{- {\sf T}} \boldsymbol{d} \Vert ^2 , \end{eqnarray*}$$

(B10)

which can be computed stably as well, by performing a forward substitution and taking the norm. Therefore, we benefit from QR decomposition even if the least squares method is not involved.

APPENDIX C: BAO MEASUREMENTS OF MOCKS WITH DIFFERENT FITTING RANGES

We perform BAO fits to the mean clustering measurements of the approximate mocks to check how fitting ranges affect BAO measurements. The corresponding covariance matrices are all rescaled by 1/1000, for better examinations on the performance of the models, as the parameters are generally well constrained with the rescaling. Thus, the priors of the parameters for this study are all flat in ranges that are sufficiently large. The bias of the α parameter is then defined as

$$\begin{eqnarray*} \delta \alpha = \frac{\alpha _{\rm fit, med} - \alpha _{\rm exp}}{ \sqrt{N_{\rm m}}\, \sigma _\alpha } , \end{eqnarray*}$$

(C1)

where N_m = 1000, α_{fit, med} indicates the median value drawn from the fitted α posterior distribution, and α_exp denotes the expected α value measured from the mocks, as is shown in Table 7.

We vary the lower and upper bounds of the fitting range in (40, 80) and |$(130, 190)\, h^{-1}\, {\rm Mpc}$|⁠, respectively, both with the step size of |$5\, h^{-1}\, {\rm Mpc}$|⁠. Using models listed in Table 6, the resulting biases and fitted errors of α, from galaxy–galaxy (ξ_gg), galaxy–void (ξ_gv), and void–void (ξ_vv) correlations of different samples, are shown in Fig. C1. In general, results for ξ_gg and ξ_gv are stable for a large set of fitting ranges, and the measured α are not significantly biased. However, BAO measurements from ξ_vv are more sensitive to the fitting ranges, especially for the lower bound (s_min). This is possibly because of the strong void exclusion patterns on small scales (see e.g. Fig. 5). Similar trends are found for the cross-correlations between the ‘LRG (c)’ and ‘ELG’ samples, as is shown in Fig. C2. Based on these results, the fiducial fitting ranges we choose are (50, 160), (60, 160), and |$(70, 160)\, h^{-1}\, {\rm Mpc}$|⁠, for ξ_gg, ξ_gv, and ξ_vv, respectively, for all the samples. In order to investigate whether the sensitivity of ξ_vv results to the fitting range bias the multitracer BAO measurements, we analyse two sets of stacked data vectors in Section 4, with or without ξ_vv included.

Figure C1.

The biases (equation C1) and fitted errors of α, from fits to the mean clustering statistics of 1000 mock realizations. In particular, for illustration purposes, the errors on α in each row are normalized with the median σ_α from the corresponding galaxy auto correlation function (ξ_gg) over all fitting ranges. Black regions indicate results outside the ranges of the colour maps, and open squares denote the fiducial fitting ranges used in this work.

Figure C2.

Similar to Fig. C1, but the results are for the cross-correlations between ‘LRG (c)’ and ‘ELG’ samples. Here, |$\xi _{\rm gg}^\times$|⁠, |$\xi _{\rm gv}^\times$|⁠, |$\xi _{\rm vg}^\times$|⁠, and |$\xi _{\rm vv}^\times$|⁠, are for the cross-correlations between LRGs and ELGs, LRGs and ELG voids, and LRG voids and ELGs, as well as LRG voids and ELG voids, respectively. The errors on α are all normalized by the median σ_α of |$\xi _{\rm gg}^\times$|⁠, over all fitting ranges. The fiducial fitting ranges are indicated by open squares.

For the multitracer approach with the combined correlation functions, we present the results with different fitting ranges and void weights in Fig. C3. The two weights are chosen to be close to the boundaries of the weight range studied in this work, i.e. (− 0.2, 0). We find that the fitting results are more sensitive to the fitting ranges when the absolute weight value gets larger, especially for the ‘LRG (c)’ sample with small s_min. This can be explained by the fact that the combined 2PCF becomes closer to the void auto correlation, and the BAO model expressed in equation (23) is not able to describe the broad-band void exclusion pattern precisely. Moreover, the void exclusion effect is more significant for ‘LRG (c)’ on the scales of interest, due to the larger void radii than the other samples. We then choose a fitting range of |$(60,150)\, h^{-1}\, {\rm Mpc}$| for all the samples, with which the fitted results are not significantly biased.

Figure C3.

Similar to Fig. C1, but the results are for the combined correlation functions evaluated using galaxy–galaxy, galaxy–void, and void–void pair counts, with two different weights applied to voids, i.e. −0.05 and −0.2. The fitted errors are normalized by the median σ_α with w = −0.05, for each sample. Open squares show our fiducial fitting ranges.

APPENDIX D: BAO MEASUREMENTS OF MOCKS WITH DIFFERENT PARAMETER PRIORS

We explore in this section impacts of parameter priors on the fitted BAO scale α, as well as its uncertainty inferred from the sampled posterior. This is done by performing BAO fits to all individual mock realizations with different parameter priors. We examine the distribution of the best-fitting (median) values of α from individual fits (α_{fit, med}), and measure the median value 〈α〉, as well as the |$1\, \sigma$| dispersion α_1σ, which is defined as half the width between the 16th and 84th percentiles of the cumulative distribution of α_{fit, med}. This dispersion is then compared to the median value of the fitted error, 〈σ_α〉, obtained from individual mocks. The relative difference between these two quantities is defined as

$$\begin{eqnarray*} \delta \sigma = \frac{\langle \sigma _{\alpha }\rangle - \alpha _{1\sigma }}{\alpha _{1\sigma }} . \end{eqnarray*}$$

(D1)

In addition, the pull statistic g(α), measured following equation (37), is also used for the prior selection.

For a given parameter p, we define its prior p_prior based on the best-fitting value (⁠|$\hat{p}$|⁠) and fitted error (⁠|$\hat{\sigma }_{p}$|⁠), obtained from fits to the mean 2PCFs of all mocks, with the covariance matrices rescaled by 1/1000. In particular, priors for fitting ξ_gg, ξ_gv, and ξ_vv of individual mocks, are all inferred from the joint fitting results with the data vector |$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}, \boldsymbol{\xi }_{\rm vv} \rbrace$|⁠, where the clustering statistics are the mean results from 1000 mock realizations. Note also that |$\hat{\sigma }_{p}$| is not rescaled here, so it is typically very small. The following priors are then tested for different parameters:

• Flat prior with a range that is considered sufficiently large;

• Flat prior |$\mathcal {F} (\hat{p} \pm n \hat{\sigma }_{p})$|⁠, where |$\mathcal {F} (p_{\rm cen} \pm p_{\rm width})$| denotes a flat prior between (p_cen − p_width) and (p_cen + p_width), and n is a fixed number;

• Flat prior |$\mathcal {U} (\hat{p} + n \hat{\sigma }_{p})$|⁠, where |$\mathcal {U} (p_{\rm max})$| indicates a flat prior between 0 and p_max;

• Gaussian prior |$\mathcal {N} (\hat{p}, n \hat{\sigma }_{p})$|⁠, where |$\mathcal {N} (p_{\rm cen}, p_{\rm sig})$| denotes a normal distribution centred at p_cen, and with the standard deviation of p_sig;

• Fixed value of |$\hat{p}$|⁠.

In particular, we fix always the c parameter in equation (23) for fitting the combined 2PCFs. The (1 + ck²) term is then equivalent to a template ‘non-wiggle’ power spectrum. By determining the value of the c parameter from the mean 2PCFs of the mocks, we essentially generate approximate templates for the combined 2PCFs with different void weights. Variu et al. (in preparation) show that this treatment enhances the robustness of the BAO measurement, without introducing significant systematic bias, compared to the result with a free c parameter. Nevertheless, we account for the potential systematic bias by varying the fixed c value (see Section 4.3).

The fitting results with different priors are listed in Tables D1 and D2, where comprehensive tests of the ‘ELG’ sample are shown, as the results are sensitive to the choices of priors, due to the large 2PCF variances of the ‘ELG’ sample. In general, the best-fitting values can be biased when no prior is imposed, and the fitted errors are significantly larger than the dispersions of best-fitting α, especially for measurements with large statistical errors (e.g. void–void clustering of the ‘ELG’ sample). This indicates that some restrictions of the parameters are necessary for improving the goodness of fits and ensuring robust error estimates. Then, we choose ranges of the priors that make the standard deviation of g(α) as close to 1 as possible. When the σ_g(α) values are similar, we choose the one that minimize the differences between the rescaled α_1σ and 〈σ_α〉.

Table D1.

Fitting results of individual mocks with different priors. Numbers in the ‘clustering’ column denote fitted errors of the mean 2PCFs of mocks, which are rescaled by |$\sqrt{1000}$|⁠. ‘–’ indicate sufficiently wide flat priors, and |$\mathcal {F}$|⁠, |$\mathcal {U}$|⁠, and |$\mathcal {N}$| denotes different prior types detailed in Appendix  D. 〈α〉 and α_1σ are the median value and |$1\, \sigma$| dispersion for the best-fitting α distribution of all mocks. 〈σ_α〉 and 〈χ²〉 denote the median value of σ_α and best-fitting χ² of the fits. δσ is the relative difference between |$\alpha _{1\, \sigma }$| and 〈σ_α〉. |$\bar{g}(\alpha)$| and σ_g(α) indicate the mean value and standard deviation of the pull quantity. Fiducial priors and the ones used for assessing systematic errors are indicated by ‘✦’ and ‘◊’, respectively.

Sample	Clustering	Bprior	Σnl, prior	〈α〉	α1σ	〈σα〉	δσ	|$\bar{g}(\alpha)$|	σg(α)	〈χ2〉/degree of freedom
MD-Patchy LRG (a)	ξgg (0.0125)	–	–	1.001	0.012	0.013	|$9.2\, \%$|	−0.04	0.97	15.4/16 = 0.96
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.000	0.012	0.012	|${\bf 3.7\, \%}$|	− 0.02	0.99	16.6/17 = 0.98	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.999	0.012	0.012	|$0.7\, \%$|	−0.03	1.02	15.5/16 = 0.97	◊
	ξgv (0.0120)	–	–	1.000	0.012	0.013	|$3.1\, \%$|	−0.00	0.97	13.4/14 = 0.95
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.000	0.012	0.012	|${\bf 0.8\, \%}$|	− 0.01	0.99	14.5/15 = 0.97	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.999	0.012	0.012	|$1.3\, \%$|	−0.01	1.01	13.4/14 = 0.96	◊
	ξvv (0.0169)	–	–	1.018	0.053	0.056	|$5.1\, \%$|	−0.11	4.42	11.2/12 = 0.94
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.000	0.017	0.017	|${\bf -1.1\, \%}$|	− 0.02	1.01	12.8/13 = 0.98	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.000	0.017	0.017	|$-2.2\, \%$|	−0.02	1.01	11.6/12 = 0.97	◊
MD-Patchy LRG (b)	ξgg (0.0115)	–	–	1.003	0.011	0.012	|$3.2\, \%$|	−0.03	0.98	15.3/16 = 0.96
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.002	0.011	0.011	|${\bf 1.0\, \%}$|	− 0.02	0.99	16.3/17 = 0.96	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.001	0.011	0.011	|$-0.9\, \%$|	−0.02	1.02	15.6/16 = 0.98	◊
	ξgv (0.0115)	–	–	1.001	0.012	0.012	|$2.0\, \%$|	−0.02	0.99	13.0/14 = 0.93
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.001	0.011	0.012	|${\bf 1.2\, \%}$|	− 0.01	0.99	14.3/15 = 0.95	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.000	0.011	0.011	|$-2.2\, \%$|	−0.02	1.02	13.1/14 = 0.94	◊
	ξvv (0.0167)	–	–	1.011	0.027	0.028	|$3.8\, \%$|	−0.34	1.83	11.2/12 = 0.94
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.000	0.016	0.016	|${\bf 3.8\, \%}$|	− 0.03	0.99	12.5/13 = 0.96	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.999	0.015	0.016	|$3.8\, \%$|	−0.05	1.02	11.5/12 = 0.96	◊
EZmock LRG (c)	ξgg (0.0152)	–	–	1.003	0.016	0.016	|$0.7\, \%$|	−0.08	0.94	15.3/16 = 0.96
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.002	0.015	0.015	|${\bf -0.4\, \%}$|	− 0.03	0.99	16.6/17 = 0.98	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.001	0.015	0.015	|$-1.0\, \%$|	−0.03	1.01	15.5/16 = 0.97	◊
	ξgv (0.0169)	–	–	1.001	0.019	0.018	|$-2.1\, \%$|	−0.00	1.03	13.1/14 = 0.94
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.000	0.017	0.017	|${\bf 0.2\, \%}$|	− 0.01	1.01	14.8/15 = 0.99	✦
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.999	0.017	0.017	|$-0.4\, \%$|	0.01	1.02	13.4/14 = 0.96	◊
	4ξvv (0.0241)	–	–	1.051	0.047	0.055	|$18.5\, \%$|	−0.10	1.70	10.8/12 = 0.90
		|$\mathcal {F}(\hat{B}\pm 50\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {F}(\hat{\Sigma }_{\rm nl}\pm 50\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.996	0.016	0.018	|$14.3\, \%$|	0.06	0.87	11.6/12 = 0.97
		|$\mathcal {F}(\hat{B}\pm 5\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	0.996	0.014	0.014	|${\bf -0.1\, \%}$|	− 0.07	1.00	13.3/13 = 1.02	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+10\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.996	0.015	0.015	|$1.5\, \%$|	−0.03	0.99	12.4/12 = 1.03	◊
EZmock ELG	6ξgg (0.0402)	–	–	1.032	0.034	0.102	|$196.9\, \%$|	−0.00	0.58	15.3/16 = 0.96
		–	|$\hat{\Sigma }_{\rm nl}$|	1.008	0.035	0.055	|$56.8\, \%$|	−0.04	0.81	16.2/17 = 0.95
		|$\mathcal {F}(\hat{B}\pm 5\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.009	0.040	0.041	|${\bf 1.9\, \%}$|	− 0.10	1.06	16.7/17 = 0.98	✦
		|$\mathcal {N}(\hat{B},3\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.009	0.040	0.041	|$1.8\, \%$|	−0.10	1.06	17.0/17 = 1.00
		|$\mathcal {F}(\hat{B}\pm 5\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {N}(\hat{\Sigma }_{\rm nl},5\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.010	0.040	0.041	|$3.8\, \%$|	−0.10	1.04	16.7/16 = 1.04
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+10\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.009	0.039	0.041	|$4.6\, \%$|	−0.10	1.03	16.1/16 = 1.01	◊
	6ξgv (0.0424)	–	–	1.027	0.037	0.099	|$163.9\, \%$|	−0.01	0.56	12.2/14 = 0.87
		–	|$\hat{\Sigma }_{\rm nl}$|	1.007	0.036	0.064	|$77.6\, \%$|	−0.04	0.79	14.1/15 = 0.94
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.008	0.044	0.044	|${\bf -1.4\, \%}$|	− 0.10	1.02	14.4/15 = 0.96	✦
		|$\mathcal {N}(\hat{B},10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.008	0.043	0.044	|$2.2\, \%$|	−0.09	1.00	14.5/15 = 0.97
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {N}(\hat{\Sigma }_{\rm nl},5\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.007	0.043	0.043	|$-0.3\, \%$|	−0.09	1.01	14.3/14 = 1.02
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+15\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.008	0.043	0.044	|$4.4\, \%$|	−0.09	1.01	13.6/14 = 0.97	◊
	6ξvv (0.0714)	–	–	1.073	0.034	0.098	|$188.1\, \%$|	0.08	0.63	10.2/12 = 0.85
		–	|$\hat{\Sigma }_{\rm nl}$|	1.027	0.045	0.089	|$98.7\, \%$|	−0.01	0.73	12.3/13 = 0.94
		|$\mathcal {F}(\hat{B}\pm 3\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.035	0.057	0.065	|${\bf 14.9\, \%}$|	− 0.04	0.91	12.7/13 = 0.98	✦
		|$\mathcal {N}(\hat{B},\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.035	0.056	0.065	|$15.0\, \%$|	−0.03	0.91	13.0/13 = 1.00
		|$\mathcal {F}(\hat{B}\pm 3\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {N}(\hat{\Sigma }_{\rm nl},5\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.033	0.054	0.065	|$19.4\, \%$|	−0.03	0.89	12.6/12 = 1.05
		|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+15\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.032	0.055	0.067	|$23.6\, \%$|	−0.06	0.98	11.6/12 = 0.96	◊
EZmock LRG (c) × ELG	2|$\xi _{\rm gg}^\times$| (0.0380)	|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.008	0.039	0.039	|${\bf -0.3\, \%}$|	− 0.08	1.07	16.2/17 = 0.96	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.008	0.037	0.039	|$5.2\, \%$|	−0.07	1.04	15.5/16 = 0.97	◊
	2|$\xi _{\rm gv}^\times$| (0.0513)	|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.008	0.051	0.055	|${\bf 8.9\, \%}$|	− 0.06	0.96	14.3/15 = 0.95	✦
		|$\mathcal {F}(\hat{B}\pm 5\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+15\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.007	0.051	0.054	|$6.1\, \%$|	−0.07	1.01	13.8/14 = 0.98	◊
	|$\xi _{\rm vg}^\times$| (0.0550)	|$\mathcal {F}(\hat{B}\pm 3\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.011	0.050	0.054	|${\bf 6.4\, \%}$|	− 0.05	0.96	14.8/15 = 0.99	✦
		|$\mathcal {F}(\hat{B}\pm 5\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+10\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.010	0.049	0.051	|$5.0\, \%$|	−0.07	1.00	13.8/14 = 0.99	◊
	|$\xi _{\rm vv}^\times$| (0.0574)	|$\mathcal {F}(\hat{B}\pm 3\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.008	0.053	0.060	|${\bf 14.2\, \%}$|	− 0.07	1.00	12.7/13 = 0.98	✦
		|$\mathcal {F}(\hat{B}\pm 5\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+15\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.009	0.052	0.062	|$18.0\, \%$|	−0.07	1.00	11.8/12 = 0.99	◊

Table D2.

Fitting results for the combined correlation functions of individual mocks, with different parameter priors. The c parameter in equation (23) is always fixed to the value for the corresponding sample in Table 6. σ_{α, mean} denotes the rescaled fitted errors with the mean 2PCFs of all mocks, and the rescaling factor is |$\sqrt{1000}$|⁠. The other notations are the same as in Table D1.

Sample	σα, mean	Bprior	Σnl, prior	〈α〉	α1σ	〈σα〉	δσ	|$\bar{g}(\alpha)$|	σg(α)	〈χ2〉/degree of
										freedom
MD-Patchy LRG (a)	0.0110	|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	0.999	0.011	0.011	|${\bf 0.7\, \%}$|	− 0.02	1.01	11.4/12 = 0.95	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.000	0.011	0.011	|$1.3\, \%$|	−0.01	1.00	12.2/13 = 0.94	◊
MD-Patchy LRG (b)	0.0107	|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+10\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.000	0.010	0.010	|${\bf 0.6\, \%}$|	− 0.02	1.04	11.6/12 = 0.97	✦
		|$\mathcal {F}(\hat{B}\pm 15\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.001	0.010	0.010	|$2.9\, \%$|	−0.02	0.99	12.5/13 = 0.96	◊
EZmock LRG (c)	0.0136	|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+30\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.002	0.014	0.014	|${\bf -1.8\, \%}$|	− 0.02	1.00	11.8/12 = 0.98	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.003	0.014	0.014	|$-3.0\, \%$|	−0.01	0.99	12.5/13 = 0.96	◊
EZmock ELG	0.0368	|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.010	0.039	0.039	|${\bf 0.6\, \%}$|	− 0.07	1.07	11.9/12 = 0.99	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.010	0.039	0.039	|$0.3\, \%$|	−0.06	0.98	12.4/13 = 0.95	◊
EZmock LRG (c) × ELG	0.0346	|$\mathcal {F}(\hat{B}\pm 10\hat{\sigma }_{_{\rm B}})$|	|$\mathcal {U}(\hat{\Sigma }_{\rm nl}+20\hat{\sigma }_{_{\Sigma _{\rm nl}}})$|	1.004	0.037	0.038	|${\bf 0.6\, \%}$|	− 0.08	1.06	11.6/12 = 0.97	✦
		|$\mathcal {F}(\hat{B}\pm 20\hat{\sigma }_{_{\rm B}})$|	|$\hat{\Sigma }_{\rm nl}$|	1.004	0.037	0.038	|$1.3\, \%$|	−0.06	1.01	12.0/13 = 0.92	◊

To be more specific, fixing the Σ_nl parameter, as is done in most BAO studies (e.g. Xu et al. 2012; Alam et al. 2017; Bautista et al. 2021; de Mattia et al. 2021), mitigates the biases of α_med in all cases. Nevertheless, without a prior of the B parameter, the median value of σ_α are still systematically larger than the dispersion of α_{fit, med}, and the distributions of g(α) do not agree with a standard normal distribution when the uncertainties are large. This is consistent with the findings in Vargas-Magaña et al. (2014). Both flat and Gaussian priors of the B parameter can solve this problem. To be more conservative, we choose flat priors as our fiducial fitting scheme. Besides, Gaussian priors of Σ_nl does not show improvements on the fitted results, but increases slightly the reduced chi-squared values. Thus, unvaried Σ_nl values are preferred in general.

When fitting to the observational data, fixed Σ_nl values require accurate estimates of the nonlinear BAO damping effect, as Σ_nl quantifies the broadening of the BAO peak. This is however only available for galaxy auto correlations ξ_gg, for which the Σ_nl value can be measured from N-body simulations. Since the BAO peak of approximate mocks is generally more damped, fitting results from the mean 2PCFs of approximate mocks can be used to estimate the upper limit of the Σ_nl value, for correlation functions involving voids. Therefore, we check also flat Σ_nl priors in the form of |$\mathcal {U} (\Sigma _{\rm nl, max})$|⁠, with only the upper limit estimated from mocks. By adjusting a bit the priors on B, the results from this type of Σ_nl priors are roughly as good as traditional BAO fitting method, i.e. with fixed Σ_nl parameters, for both individual correlation functions and the combined 2PCF.

We then further check the BAO measurements with multiple clustering measurements fitted simultaneously, and the results are listed in Table D3. This time, the pull quantities with fixed Σ_nl prior are still consistent with those of the standard normal distribution. However, for the Σ_nl prior with an upper limit, the fitted errors appear to be underestimated for the ‘ELG’ sample, as well as the cross-correlation measurements between the ‘LRG (c)’ and ‘ELG’ samples. This is possibly because the clustering measurements in these cases are too noisy to have constraints on Σ_nl. Thus, eventually we use fixed Σ_nl values for fits to multiple 2PCFs. Nevertheless, accurate estimations of Σ_nl are crucial for the combined correlation functions, as there is no way so far to measure the values using N-body simulations. Since the Σ_nl priors with upper limits work well for the combined 2PCFs (see Table D2), we choose these priors as our fiducial fitting scheme for the combined correlation functions. Moreover, we compare the results with the two types of Σ_nl priors, to assess the potential systematic biases (see Section 4.3).

Table D3.

Multitracer fitting results for combinations of correlation functions of individual mocks, with different parameter priors. Rows marked with ‘✦’ and ‘◊’ are obtained with priors indicated by the same symbol in Table D1. Meanings of columns with numbers are the same as in Tables D1 and D2.

Sample	Data vector	σα, mean	Bprior	Σnl, prior	〈α〉	α1σ	〈σα〉	δσ	|$\bar{g}(\alpha)$|	σg(α)	〈χ2〉/degree of
											freedom
MD-Patchy LRG (a)	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}\rbrace$|	0.0112	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.001	0.011	0.011	|${\bf 3.9\, \%}$|	− 0.01	0.98	32.1/33 = 0.97	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.000	0.011	0.010	|$-5.0\, \%$|	−0.01	1.05	29.9/31 = 0.96	◊
	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}, \boldsymbol{\xi }_{\rm vv}\rbrace$|	0.0112	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.001	0.011	0.011	|${\bf 0.5\, \%}$|	− 0.01	1.00	45.0/47 = 0.96	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.000	0.011	0.010	|$-8.9\, \%$|	−0.01	1.08	42.3/44 = 0.96	◊
MD-Patchy LRG (b)	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}\rbrace$|	0.0104	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.002	0.010	0.010	|${\bf 4.0\, \%}$|	− 0.02	0.99	31.9/33 = 0.97	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.002	0.010	0.010	|$-2.9\, \%$|	−0.02	1.07	30.3/31 = 0.98	◊
	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}, \boldsymbol{\xi }_{\rm vv}\rbrace$|	0.0106	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.002	0.010	0.010	|${\bf 0.6\, \%}$|	− 0.02	1.02	45.0/47 = 0.96	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.001	0.010	0.009	|$-6.9\, \%$|	−0.02	1.10	43.0/44 = 0.98	◊
EZmock LRG (c)	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}\rbrace$|	0.0144	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.002	0.014	0.014	|${\bf -1.5\, \%}$|	− 0.02	1.02	32.6/33 = 0.99	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.002	0.015	0.014	|$-7.2\, \%$|	−0.01	1.08	30.2/31 = 0.97	◊
	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}, \boldsymbol{\xi }_{\rm vv}\rbrace$|	0.0136	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.001	0.011	0.011	|${\bf -2.0\, \%}$|	− 0.04	1.02	46.3/47 = 0.99	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	0.999	0.011	0.011	|$-4.9\, \%$|	−0.02	1.07	42.9/44 = 0.97	◊
EZmock ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}\rbrace$|	0.0376	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.012	0.039	0.038	|${\bf -3.3\, \%}$|	− 0.09	1.02	31.8/33 = 0.96	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.010	0.040	0.037	|$-7.0\, \%$|	−0.09	1.09	30.7/31 = 0.99	◊
	|$\lbrace \boldsymbol{\xi }_{\rm gg}, \boldsymbol{\xi }_{\rm gv}, \boldsymbol{\xi }_{\rm vv}\rbrace$|	0.0379	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.014	0.040	0.037	|${\bf -6.9\, \%}$|	− 0.11	1.09	45.3/47 = 0.96	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.012	0.042	0.035	|$-17.9\, \%$|	−0.12	1.26	42.8/44 = 0.97	◊
EZmock LRG (c)× ELG	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times , \boldsymbol{\xi }_{\rm gv}^\times , \boldsymbol{\xi }_{\rm vg}^\times \rbrace$|	0.0370	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.008	0.037	0.036	|${\bf -2.1\, \%}$|	− 0.08	1.06	46.3/49 = 0.95	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.007	0.039	0.031	|$-18.7\, \%$|	−0.10	1.20	44.0/46 = 0.96	◊
	|$\lbrace \boldsymbol{\xi }_{\rm gg}^\times , \boldsymbol{\xi }_{\rm gv}^\times , \boldsymbol{\xi }_{\rm vg}^\times , \boldsymbol{\xi }_{\rm vv}^\times \rbrace$|	0.0370	|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \hat{\Sigma }_{\rm nl} \rbrace$|	1.009	0.036	0.035	|${\bf -2.3\, \%}$|	− 0.10	1.09	59.0/63 = 0.94	✦
			|$\lbrace \mathcal {F}\rbrace$|	|$\lbrace \mathcal {U}\rbrace$|	1.009	0.038	0.030	|$-21.6\, \%$|	−0.15	1.32	56.1/59 = 0.95	◊

With the fiducial fitting priors, the median values of α, as well as their dispersions, are generally consistent with the results from fits to the mean 2PCFs of all mock realizations. Obvious discrepancies are only observed for few cases, such as for the void auto correlations of the ‘LRG (c)’ sample. In this case, the pull distribution, as well as the comparison of the dispersion of α_{fit, med} and the fitted σ_α, both prefer the smaller BAO uncertainty from fits to individual mocks. Even with prior widths of 100 times |$\hat{\sigma }_{_B}$| and |$\hat{\sigma }_{_{\Sigma _{\rm nl}}}$|⁠, the fitted σ_α from individual mocks are still significantly smaller than that of the mean 2PCFs of all mocks. Besides, we have further performed a fit to the mean ξ_vv of all mocks, with our fiducial fitting priors, and the result is consistent with that of individual mocks. It implies that a prior is always necessary for void clustering of the ‘LRG (c)’ sample. There is however no such problem for the combined 2PCFs. This suggests that a multitracer BAO fitting approach based on the combined 2PCFs is more robust, than that relying on the stacked data vectors of multiple 2PCFs.

APPENDIX E: BAO MEASUREMENTS FOR DIFFERENT GALACTIC CAPS

The 2PCFs of SDSS samples in different Galactic caps, as well as the 1σ regions of the corresponding approximate mocks, are shown in Fig. E1. It can be seen that measurements from both Galactic caps are generally consistent, though the differences for those involving ELGs are relatively large (see also Raichoor et al. 2021). To examine the statistical consistencies of BAO measurements from different Galactic caps, we analyse the SDSS data in NGC and SGC separately with the fitting scheme shown in Table 6. The resulting marginalized posterior distributions of α are shown in Fig. E2, in comparison with results from the combined NGC and SGC samples. Here, the galaxy-only results for the eBOSS samples except ELGs in NGC are consistent with those from the joint BAO and RSD analysis in Wang et al. (2020). However, the posterior distributions of α measured from the ‘ELG’ sample in NGC seem to prefer values outside our prior range, for both galaxy auto correlation and the combined 2PCF of galaxies and voids. This can be explained by the failure of BAO detection from this subsample, which is consistent with the relatively large cosmic variance of the ‘ELG’ sample, based on the analyses of individual mocks (Raichoor et al. 2021).

Figure E1.

Galaxy–galaxy and combined 2PCFs measured from samples in different Galactic caps. Envelopes show the 1σ dispersions of 2PCFs from mocks. The abscissas of the SDSS NGC and SGC measurements are shifted by −0.5 and |$0.5\, h^{-1}\, {\rm Mpc}$|⁠, respectively for illustration purposes.

Figure E2.

Marginalized posterior distributions of α measured from the 2PCFs of SDSS data in different Galactic caps. The tensions between NGC and SGC measurements, as well as the fractions of mocks with larger tensions, are given in the subplots except for the ‘ELG’ sample.

We define the tension between BAO measurements from NGC and SGC as

$$\begin{eqnarray*} \delta \alpha _{_{\rm NS}} = \frac{| \alpha _{_{\rm NGC}} - \alpha _{_{\rm SGC}} |}{\sqrt{\sigma _{\alpha , \scriptscriptstyle {\rm NGC}}^2 + \sigma _{\alpha , \scriptscriptstyle {\rm SGC}}^2}} , \end{eqnarray*}$$

(E1)

where |$\alpha _{_{\rm NGC/SGC}}$| and |$\sigma _{\alpha , \scriptscriptstyle {\rm NGC/SGC}}$| denote the median value and 1σ errors of α respectively, which are obtained from the corresponding posterior distributions. The tensions measured from different samples are shown in Fig. E2 as well, together with the fractions of mock realizations that present larger differences between the two Galactic caps. We observe a relatively large tension (2.5σ) for the galaxy auto correlation of the ‘LRG (a)’ sample. Nevertheless, the tension is significantly reduced when including the clustering of voids, in which case almost 10 per cent of the mocks show larger tensions. For the rest of the samples, apart from ‘ELG’, the BAO measurements from the combined correlation functions in NGC and SGC are well consistent, with the differences smaller than 1σ. The results suggest that it is reasonable to combine the clustering statistics in NGC and SGC for the multitracer BAO measurements.

APPENDIX F: CROSS COVARIANCES OF DISTANCE MEASUREMENTS FROM DIFFERENT SAMPLES

Since some of our galaxy samples share the same volume, or even contain identical objects, it is crucial to take into account their cross-correlations for cosmological constraints (see Fig. A2). However, the cross covariances of distance measurements are not directly available from the BAO fits. Thus, we use the correlation matrices R of |$D_{_{\rm V}} / r_{\rm d}$| obtained from individual mock realizations to approximate the covariances of measurements from different samples. Denoting the statistical and systematic errors of distances measured from the SDSS data as |$\boldsymbol{\sigma }$| and |$\boldsymbol{s}$| respectively, where the vector components are results from individual samples, the covariance matrix for the measurements of all samples is computed as

$$\begin{eqnarray*} {\bf C} = m_1 (\boldsymbol{\sigma } \otimes \boldsymbol{\sigma }) \circ {\bf R} + \boldsymbol{s}^2 {\bf I} . \end{eqnarray*}$$

(F1)

Here, ⊗ and ○ are the outer product and Hadamard (element-wise) product respectively, and I indicates the identity matrix. The factor m₁ given by equation (32) is for accounting the precision of the covariance matrices used for BAO fits. Since we always use the same number of mocks for covariance matrix estimates, as well as the same number of parameters and data points for the fits, m₁ is identical for all samples. Furthermore, we assume no cross-correlations between the systematic errors of different samples, so they are only added to the diagonal terms of the covariance matrix.

Using the correlation coefficients shown in Fig. 14, the full covariance matrix of |$D_{_{\rm V}} / r_{\rm d}$| measured from galaxy–galaxy correlations of the five SDSS samples is

$$\begin{eqnarray*} &&{\bf C}_{\rm gg} (D_{_{\rm V}}/r_{\rm d}) =\nonumber \\ &&\begin{array}{ccccccc}\scriptstyle z_{\rm eff}: & {\scriptstyle 0.38} & {\scriptstyle 0.51} & {\scriptstyle 0.70} & {\scriptstyle 0.77} & {\scriptstyle 0.85} & \end{array}\nonumber \\ &&10^{-3} \begin{pmatrix}& 9.9313 & 4.6559 & & & &\\ & 4.6559 & 13.651 & & & & \\ & & & 37.753 & 13.863 & 6.2144 & \\ & & & 13.863 & 195.24 & 86.093 & \\ & & & 6.2144 & 86.093 & 372.50 & \\ \end{pmatrix} \begin{array}{c}\scriptstyle {\rm LRG\, (a)} \\ {\scriptstyle {\rm LRG\, (b)}}\\ {\scriptstyle {\rm LRG\, (c)}}\\ {\scriptstyle {\rm CROSS}}\\ {\scriptstyle {\rm ELG}} \end{array}, \end{eqnarray*}$$

(F2)

where the effective redshifts of different samples are shown, and ‘CROSS’ denotes the results from cross-correlations between the ‘LRG (c)’ and ‘ELG’ samples. Similarly, the covariance matrix for distances measured from the combined correlations is

$$\begin{eqnarray*} &&{\bf C}_{\rm comb} (D_{_{\rm V}}/r_{\rm d}) =\nonumber \\ &&\begin{array}{ccccccc}\scriptstyle z_{\rm eff}: & {\scriptstyle 0.38} & {\scriptstyle 0.51} & {\scriptstyle 0.70} & {\scriptstyle 0.77} & {\scriptstyle 0.85} & \\ \end{array}\nonumber \\ &&10^{-3}\begin{pmatrix}& 7.2216 & 3.4095 & & & & \\ & 3.4095 & 9.8147 & & & & \\ & & & 33.825 & 12.042 & 3.4234 & \\ & & & 12.042 & 149.02 & 61.581 & \\ & & & 3.4234 & 61.581 & 278.53 & \\ \end{pmatrix} \begin{array}{l}\scriptstyle {\rm LRG\, (a)}\\ {\scriptstyle {\rm LRG\, (b)}}\\ {\scriptstyle {\rm LRG\, (c)}}\\ {\scriptstyle {\rm CROSS}}\\ {\scriptstyle {\rm ELG}}\\ \end{array} . \end{eqnarray*}$$

(F3)