Cosmological Constraints from the BOSS DR12 Void Size Function

The comoving number counts of voids as a function of their radius, i.e., the void size function, has been modeled from first principles by Sheth & van de Weygaert (2004). This model is developed within the framework of the excursion-set theory, with the same formalism adopted for the halo mass function (Press & Schechter 1974; Peacock & Heavens 1990; Bond et al. 1991; Cole 1991; Mo & White 1996). The multiplicity function, expressing in this case the volume fraction of the Universe occupied by cosmic voids in linear theory, ⁸ yields

$$\begin{eqnarray}&&{f}_{\mathrm{ln}{\sigma }_{{\rm{m}}}}({\sigma }_{{\rm{m}}})=2\displaystyle \sum _{j=1}^{\infty }\,\exp \left(-\displaystyle \frac{{\left(j\pi x\right)}^{2}}{2}\right)\,j\pi {x}^{2}\,\sin \left(j\pi { \mathcal D }\right),\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}&&{ \mathcal D }=\displaystyle \frac{| {\delta }_{{\rm{v}}}^{{\rm{L}}}| }{{\delta }_{{\rm{c}}}^{{\rm{L}}}+| {\delta }_{{\rm{v}}}^{{\rm{L}}}| },\qquad x=\displaystyle \frac{{ \mathcal D }}{| {\delta }_{{\rm{v}}}^{{\rm{L}}}| }{\sigma }_{{\rm{m}}},\end{eqnarray} \tag{ 2 }$$

where σ_m is the rms variance of the linear matter perturbations on a scale R_L, while ${\delta }_{{\rm{v}}}^{{\rm{L}}}$ and ${\delta }_{{\rm{c}}}^{{\rm{L}}}$ are the two fundamental barriers of the model. ⁹

The multiplicity function of Equation (1) assumes initial spherical density fluctuations to evolve according to linear theory over cosmic time. Therefore, to predict the number counts of observed voids of size R, we need to consider the nonlinear evolved density field. For this purpose, Jennings et al. (2013) modified the void size function model imposing the conservation of the volume occupied by voids during the transition from linearity to nonlinearity, which for this reason is called volume-conserving (hereafter, Vdn) model and is expressed as

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dn}}{d\mathrm{ln}R}=\displaystyle \frac{{f}_{\mathrm{ln}{\sigma }_{{\rm{m}}}}({\sigma }_{{\rm{m}}})}{V(R)}\,\displaystyle \frac{d\mathrm{ln}{\sigma }_{{\rm{m}}}^{-1}}{d\mathrm{ln}{R}_{{\rm{L}}}}\right|}_{{R}_{{\rm{L}}}={R}_{{\rm{L}}}(R)}.\end{eqnarray} \tag{ 3 }$$

The Vdn model successfully predicts the size function of cosmic voids identified in the total matter field (i.e., in the distribution of dark matter particles), provided that these voids are defined as spherical nonoverlapping spheres embedding a specific density contrast ${\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}}$ (Jennings et al. 2013; Ronconi & Marulli 2017; Contarini et al. 2019; Ronconi et al. 2019). The threshold ${\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}}$ is the nonlinear counterpart of the quantity ${\delta }_{{\rm{v}}}^{{\rm{L}}}$ introduced in Equation (2), where the subscript "DM" specifies the requirement to be computed in the dark matter density field. To convert the nonlinear threshold into its linear counterpart and vice versa, we can make use of the fitting formula provided by Bernardeau (1994):

$$\begin{eqnarray}&&{\delta }_{{\rm{v}}}^{{\rm{L}}}={ \mathcal C }\,\left[1-{\left(1+{\delta }_{{\rm{v}}}^{\mathrm{NL}}\right)}^{-1/{ \mathcal C }}\right],\mathrm{with}\ { \mathcal C }=1.594,\end{eqnarray} \tag{ 4 }$$

which features enhanced precision when applied to underdensities.

One of the observational effects to consider when computing the theoretical void size function is the so-called geometric distortions, which are caused by the assumption of a fiducial cosmology that may not coincide with the true one. Geometric distortions affect the observed void shape by introducing an anisotropy between the directions parallel and perpendicular to the line of sight. This distortion is commonly denoted as the Alcock & Paczynski (1979; AP) effect and can be modeled following the prescriptions by, e.g., Sánchez et al. (2017a):

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\parallel } & = & \displaystyle \frac{{H}^{* }(z)}{H(z)}\,{r}_{\parallel }^{* }\equiv {q}_{\parallel }\,{r}_{\parallel }^{* },\\ {r}_{\perp } & = & \displaystyle \frac{{D}_{{\rm{A}}}(z)}{{D}_{{\rm{A}}}^{* }(z)}\,{r}_{\perp }^{* }\equiv {q}_{\perp }\,{r}_{\perp }^{* },\end{array}\end{eqnarray} \tag{ 5 }$$

where ${r}_{\parallel }^{* }$ and ${r}_{\perp }^{* }$ are the observed comoving distances between two objects at redshift z, projected along the parallel and perpendicular directions, respectively, H(z) is the Hubble parameter, and D_A(z) is the comoving angular-diameter distance. The quantities marked with asterisks are those computed by assuming the fiducial cosmology, while the remaining ones are calculated in the true cosmology. In the case of cosmic voids, we can derive the radius R of a sphere having the same volume of the observed deformed void as $R={q}_{\parallel }^{1/3}{q}_{\perp }^{2/3}{R}^{* }$ (Ballinger et al. 1996; Eisenstein et al. 2005; Xu et al. 2013; Sánchez et al. 2017a; Hamaus et al. 2020; Correa et al. 2021). The application of this correction to the void size function model translates into a shift of the predicted void counts toward smaller or larger radii, depending on the parameter values assumed for the fiducial cosmology.

When dealing with real survey data, cosmic voids cannot be defined in the dark matter particle field, so an extension of this modeling to biased tracers, such as galaxies and galaxy clusters, is required. In particular, it is fundamental to consider the effect of tracer bias on the void density profiles and consequently on the associated void radius.

The density profiles of voids traced by biased objects are deeper than those traced by dark matter particles. The radius corresponding to a fixed internal spherical density contrast is consequently larger for the former (e.g., Sutter et al. 2014a; Pollina et al. 2017; Schuster et al. 2023). To include this effect in the void size function theory, Contarini et al. (2019) proposed a parameterization of the Vdn model nonlinear underdensity threshold as a function of the tracer effective bias, b_eff:

$$\begin{eqnarray}\begin{array}{rcl}{\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}} & = & \displaystyle \frac{{\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}}{{ \mathcal F }({b}_{\mathrm{eff}},z)},\mathrm{with}\\ { \mathcal F }({b}_{\mathrm{eff}},z) & = & {B}_{\mathrm{slope}}\,{b}_{\mathrm{eff}}(z)+{B}_{\mathrm{offset}},\end{array}\end{eqnarray} \tag{ 6 }$$

where ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}$ is the tracer density contrast embedded by voids identified in the tracer field and, as already mentioned, ${\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}}$ is the matter underdensity threshold to be used in the void size function model, after its conversion from linear theory. The two coefficients regulating the slope and the offset of the linear function ${ \mathcal F }({b}_{\mathrm{eff}})$, i.e., B_slope and B_offset, are redshift independent and can be derived by calibrating the void size function model with simulations (Contarini et al. 2019, 2021, 2022). They have been shown to depend on the kind of tracers used to identify the voids, but also to have a negligible dependence on the cosmological model (Contarini et al. 2021).

Despite its simplicity and demonstrated effectiveness, this model parameterization requires an accurate measurement of the tracer effective bias, b_eff. This quantity is not trivial to estimate in real data catalogs, since it exhibits a degeneracy with the matter power spectrum normalization, σ₈. For example, by modeling the two-point correlation function (2PCF), only the product of these parameters, b_eff σ₈, can be derived. Therefore, the assumption of a fiducial value of σ₈ is necessary to recover separately the galaxy effective bias. To avoid this problem, in this work we introduce a different parameterization of the linear function ${ \mathcal F }$, which also depends on σ₈:

$$\begin{eqnarray}&&{ \mathcal F }({b}_{\mathrm{eff}}\,{\sigma }_{8},z)={C}_{\mathrm{slope}}\,{b}_{\mathrm{eff}}(z)\,{\sigma }_{8}(z)+{C}_{\mathrm{offset}}.\end{eqnarray} \tag{ 7 }$$

This expression has the advantage to be based on a σ₈-dependent bias, i.e., the product b_eff σ₈, which can be estimated directly via the modeling of the 2PCF.

The linear parameterization reported in Equation (6) and equivalently the form of Equation (7) is effective also to encapsulate redshift-space distortions (Contarini et al. 2022), i.e., the deformation of void shapes due to the peculiar velocities of tracers. Indeed, when estimating comoving distances from redshifts, the coherent outflow of tracers around the void center leads to larger observed void radii on average (Pisani et al. 2015b).

Therefore, in this work we model the void size function by means of the Vdn theory modified to correct for both geometric and dynamic effects, relying on an underdensity threshold parameterization dependent on the quantity b_eff σ₈. Hereafter, we will refer to this model as extended Vdn.

In this work we analyze the BOSS DR12 galaxy sample, covering both the northern and the southern galactic hemispheres, for a total sky area of more than 10,000 square degrees. We account for the two target selections denoted LOWZ and CMASS (Reid et al. 2016), whose combination features more than 1,000,000 galaxies spanning the redshift range 0.2 < z < 0.75 and constitutes one of largest spectroscopic catalogs available to date.

We also make use of the MultiDark PATCHY mock catalogs, specifically designed to mimic relevant observational effects characterizing the BOSS DR12 galaxy sample, including selection functions and masking. In particular, we consider 100 different realizations of these mocks to calibrate and validate our void size function model. These catalogs have been produced by populating the dark matter halos of the MultiDark cosmological simulations (Klypin et al. 2016) with the PATCHY method (Kitaura et al. 2014), which accurately reproduces the clustering properties of different populations of objects (Kitaura et al. 2016b). In particular, Kitaura et al. (2016a) and Rodríguez-Torres et al. (2016) demonstrated that the PATCHY mocks follow accurately the structure formation down to scales of a few megaparsecs, as well as the evolution of galaxy bias and nonlinear redshift-space distortions. Moreover, PATCHY mocks reproduce the BOSS DR12 power spectrum, two- and three-point correlation functions up to k ∼ 0.3 h Mpc⁻¹, for arbitrary stellar mass bins.

The MultiDark simulation has been run with a flat ΛCDM cosmology characterized by Planck2013 parameters (Planck Collaboration et al. 2014): Ω_m = 0.307115, Ω_b = 0.048206, σ₈ = 0.8288, n_s = 0.9611, and h = 0.6777. Although many mock realizations are available, all of the PATCHY simulations are built with this set of cosmological parameters. We do not expect this to impact significantly our results, despite our aim at analyzing cosmological models with different values of Ω_m, σ₈, and w. Indeed, the effectiveness of the void size function model has been proved across various cosmological parameter ranges and alternative cosmological scenarios (see, e.g., Verza et al. 2019; Contarini et al. 2021).

We identify the cosmic underdensities in the distribution of galaxies both in real and simulated data by means of the public Void IDentification and Examination toolkit ¹⁰ (VIDE;Sutter et al. 2015). VIDE is a parameter-free void finding algorithm based on the ZOnes Bordering On Voidness (ZOBOV; Neyrinck 2008) code, which performs a Voronoi tessellation of tracer particles to reconstruct a continuous density field. Then it finds the local density minima, from which the Voronoi cells start to be merged progressively to create density basins through a watershed algorithm (Platen et al. 2007). The whole 3D space is thereby divided up by thin overdense ridges into separate underdense regions, forming the final catalog of voids. For each underdensity identified, the corresponding void center position X is calculated as the volume-weighted barycenter of the N Voronoi cells that define the void, while the void effective radius R_eff is taken as the radius of a sphere with the total volume of all of those cells:

$$\begin{eqnarray}&&{\boldsymbol{X}}\equiv \displaystyle \frac{{\sum }_{i=1}^{N}{{\boldsymbol{x}}}_{i}{V}_{i}}{{\sum }_{i=1}^{N}{V}_{i}}\quad {\rm{}}\,{\rm{and}}\quad {R}_{\mathrm{eff}}\equiv {\left(\displaystyle \frac{3}{4\pi }\displaystyle \sum _{i=1}^{N}{V}_{i}\right)}^{1/3}{\rm{}}\,,\end{eqnarray} \tag{ 8 }$$

where x _i are the coordinates of the ith tracer of a given void, and V_i is the volume of its associated Voronoi cell.

We apply VIDE to the BOSS DR12 galaxy sample and to the 100 PATCHY realizations, assuming throughout the analysis the Planck2013 cosmological parameters. This cosmology is used to convert angles and redshifts into comoving coordinates to build the galaxy map, and is assumed as fiducial to compute the void size function model in the Bayesian approach we present in Section 3.3. The impact of the fiducial cosmology assumption on the results presented in this paper is analyzed in Appendix C.

The number of cosmic voids identified by the algorithm is strictly related to both the survey volume and the number density of the objects used to trace voids (see, e.g., Schuster et al. 2023). We report in the left panel of Figure 1 the number density of galaxies computed in concentric redshift shells, n_gal(z), for both BOSS and PATCHY samples. We show the results for the PATCHY data set as a shaded band representing the 68% and 95% confidence levels around the mean number density value, computed from the standard deviation between the 100 mock realizations. The PATCHY data standard deviation is used also to calculate the uncertainty associated with the number density of BOSS galaxies, reported in this plot with 2σ error bars for visibility, and interpolated with a spline function. As a comparison, we also show the number density distribution of mock galaxies processed by means of the weights specifically designed for the MultiDark PATCHY catalogs, which mimic observational systematic effects that affect the completeness of the galaxy sample. We underline, however, that the weighted PATCHY samples will not be used in the following analysis, since our aim is to exploit mock catalogs to calibrate the void size function model and validate our pipeline. Indeed, using galaxy catalogs characterized by a higher completeness allows us to obtain tighter constraints on the calibration parameters and does not bias the final results with BOSS data, since the galaxy clustering of the two samples is fully compatible on the relevant range of scales (see Appendix A).

Zoom In Zoom Out Reset image size

As expected, the intrinsic higher mean density of the PATCHY galaxy catalog with respect to the BOSS one is reflected in the larger void numbers of these samples. The total number of voids identified with the void finder VIDE above a minimum effective radius of 30 h⁻¹ Mpc is N_v = 3280 for the BOSS DR12 data and ${\overline{N}}_{{\rm{v}}}=5111.85$ as a mean value for the 100 PATCHY mock realizations. We note that in this estimate, as well as in the analysis we will present, we exclude galaxies with z > 0.65, since their number density is very low in this regime and extended objects like voids are potentially affected by the survey boundary.

We report in the rightmost panels of Figure 1 the size distribution of void number counts, dividing the sample in two bins of redshift, 0.2 < z ≤ 0.45 and 0.45 < z < 0.65, approximately coinciding with the two target selections LOWZ and CMASS. We show the void count size distribution from PATCHY via two shaded bands representing confidence levels of 68% and 95% as obtained from the 100 realizations. The corresponding uncertainty on the BOSS void number counts is computed by rescaling the PATCHY 68% errors by the square root of the ratio between BOSS and PATCHY survey volumes. Moreover, we use this survey volume ratio between PATCHY mocks and BOSS data to rescale the void counts in PATCHY, as shown in the two rightmost panels of Figure 1.

We computed the exact volume for these galaxy catalogs by summing up all Voronoi cells created by VIDE to tessellate the 3D space. This corresponds approximately to the total volume occupied by all of the voids identified by the algorithm. This approach also takes into account the removal of voids intersecting survey boundaries (see Sutter et al. 2015 for further details). The effective volumes derived with this procedure will be used in the following analysis also to convert the void number density predictions provided by the extended Vdn theory into void counts, N_voids, to model accordingly the measured void size function.

As is evident from Figure 1, the size distribution of the rescaled PATCHY void counts matches the BOSS data closely. Only toward the small-scale cut, at R_eff < 35 h⁻¹ Mpc, do small deviations occur, which may be caused by finite resolution effects close to the mean galaxy separation $\mathrm{MGS}{(z)={n}_{\mathrm{gal}}(z)}^{-1/3}$ (see, e.g., Jennings et al. 2013; Pisani et al. 2015a; Ronconi & Marulli 2017; Contarini et al. 2019, 2022, for further details). Following our cleaning procedure, as described in Section 3.2, these small voids will be excluded in our analysis anyway. Finally, we note that the void counts from PATCHY can reach numbers below unity due to the averaging procedure over all mock catalogs.

The void sample provided by VIDE is not yet directly exploitable for the study of the void size function (Jennings et al. 2013; Ronconi & Marulli 2017), unless a number of ad hoc correction factors are introduced in order to link the cosmic voids found by the void finder with the ideal voids considered in the void size function theory. As introduced in Section 2, the Vdn model predicts the number counts of voids defined as spherical nonoverlapping regions embedding a specific density contrast value in the total matter field, ${\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}}$.

In this work, we adopt the complementary approach of aligning the considered sample of voids with the prescriptions of the Vdn model (as previously done by Contarini et al. 2019, 2021, 2022; Ronconi et al. 2019; Verza et al. 2019), therefore without modifying the void definition used in the theoretical model. To this end, we follow the cleaning methodology proposed by Jennings et al. (2013), applying an upgraded version of the code implemented by Ronconi & Marulli (2017). In this algorithm, the only settable parameter is the threshold ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}$, which is computed in the density field of tracers (i.e., galaxies in this case) and is used to prepare the void sample. In this upgraded version of the cleaning algorithm, the galaxy density distribution is approximated as reported in the left panel of Figure 1, i.e., by a spline fit to the number density of galaxies as a function of redshift.

In our analysis, we follow the choice of Contarini et al. (2022) to select a threshold of ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}=-0.7$: this density contrast has been proven to represent a good compromise to model the innermost underdense regions of voids and simultaneously avoid poor statistics due to the sparsity of tracers. We underline that the selection of this threshold is not unique. We refer the reader to Contarini et al. (2022) for further details about the impact of this parameter choice.

Once the threshold ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}$ has been selected, the cleaning algorithm first discards from the original VIDE void catalog all of those underdensities having an internal spherical density contrast excessively high: starting from the void center identified by VIDE, ¹¹ a sphere with starting radius R_central = 2 MGS(z) is grown, and if this sphere does not match the minimum requirement imposed by the chosen threshold, i.e., ${\delta }_{{\rm{v}},\mathrm{tr}}({R}_{\mathrm{central}})\lt -0.7$, the corresponding void is rejected. Second, all of the remaining underdensities are rescaled following the shape of their void density profile up to a radius R, ¹² embedding a spherical density contrast equal to ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}$. Lastly, all of the mutually intersecting voids are removed by rejecting, one by one, the overlapping object with the higher central density.

The final outcome of this cleaning procedure is a sample of voids ideally free of spurious Poissonian fluctuations and shaped in agreement with the theoretical definition adopted in the void size function model. In particular, to match the void number counts measured from this sample and the void size function predicted by the extended Vdn model (see Section 2), we must convert the underdensity threshold used to rescale the voids, ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}=-0.7$ in our case, by means of Equation (7); this allows us to recover the corresponding value in the total matter density field, ${\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}}$, which then is inserted in Equation (2) after the transformation to its linear theory counterpart via Equation (4).

We show an example of the cleaned void samples obtained in this work in Figure 2, where we report 200 h⁻¹ Mpc-thick slices of the BOSS DR12 and MultiDark PATCHY light-cones. In this visual representation, we depict galaxies and voids in the full redshift range 0.2 < z < 0.75, separating the observed BOSS data (top-left panels) from those corresponding to one of the 100 PATCHY mocks (bottom-right panels). As expected, we see a similarity between the projected distribution of galaxies and voids in the specular sky sections of the observed and simulated data. We also point out that any apparent overlapping between voids is given by projection effects within the slice and that the regions not featuring the presence of voids are those excluded by the survey mask.

Zoom In Zoom Out Reset image size

Now that the void size function theory and the void sample preparation have been outlined, we describe the Bayesian analysis we perform to investigate different cosmological scenarios and to calibrate the void size function model (see Section 3.4). According to Bayesian statistics, the posterior probability of a set of model parameters, Θ, given a data set, ${ \mathcal D }$, which is composed of the measured void number counts in our case, can be computed as ${ \mathcal P }({\rm{\Theta }}| { \mathcal D })\propto { \mathcal L }({ \mathcal D }| {\rm{\Theta }})\,p({\rm{\Theta }})$, where ${ \mathcal L }({ \mathcal D }| {\rm{\Theta }})$ is the likelihood, and p(Θ) is the prior distribution. Since in this work we consider the number counts of cosmic voids, the likelihood can be assumed to follow Poisson statistics (Sahlén et al. 2016):

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }({ \mathcal D }| {\rm{\Theta }})\\ \quad =\,\displaystyle \prod _{i,j}\displaystyle \frac{N{\left({r}_{i},{z}_{j}| {\rm{\Theta }}\right)}^{N({r}_{i},{z}_{j}| { \mathcal D })}\,\exp \left[-N({r}_{i},{z}_{j}| {\rm{\Theta }})\right]}{N({r}_{i},{z}_{j}| { \mathcal D })!},\end{array}\end{eqnarray} \tag{ 9 }$$

where the product is over the radius and redshift bins, labeled with i and j, respectively. The quantity $N({r}_{i},{z}_{j}| { \mathcal D })$ is the number of voids in the ith radius bin and jth redshift bin, while N(r_i , z_j ∣Θ) corresponds to the expected value in the considered cosmological model, given a set of parameters Θ.

The first cosmological scenario we investigate in this work is the standard flat ΛCDM model, characterized by six fundamental parameters: the total matter density parameter (Ω_m), the dimensionless Hubble parameter (h ≡ H₀/100 km s⁻¹ Mpc⁻¹), the baryon density parameter times h² (Ω_b h²), the spectral index of the primordial power spectrum (n_s), the reionization optical depth (τ), and the primordial power spectrum amplitude (A_s), which we appropriately reinterpret with the present-time matter power spectrum normalization (σ₈). In our modeling, we require flatness by imposing the cosmological constant-density parameter to be Ω_Λ = 1 − Ω_m and define the value of the cold dark matter component as Ω_cdm = Ω_m − Ω_b. We underline that, analogously to σ₈, the Universe density parameters defining the cosmological scenario are expressed with their present-day values.

In addition to the current standard model of cosmology, we also consider a Universe framework implementing a constant dark energy equation of state, which is commonly denoted wCDM. This model represents the simplest dark energy model alternative to ΛCDM, and is therefore characterized by an additional parameter, w. This scenario reduces to standard ΛCDM when w = −1.

In the Bayesian analysis reported in Section 4, both cosmological scenarios will be explored. We focus in particular on the cosmological degeneracies acting on the parameter planes Ω_m–σ₈, for the ΛCDM model, and Ω_m–w, for the wCDM model. This selection is motivated by recent promising studies concerning the complementarity of different void statistics with standard cosmological probes, which showed the effectiveness of a joint analysis in constraining the growth of cosmic structures and the dark energy equation of state (Pisani et al. 2015a; Bayer et al. 2021; Bonici et al. 2022; Contarini et al. 2022; Hamaus et al. 2022; Kreisch et al. 2022; Pelliciari et al. 2023).

We sample the posterior distribution of the model parameters with a Markov Chain Monte Carlo (MCMC) technique, by considering wide uniform priors for the cosmological parameters of interest, i.e., Ω_m and σ₈ for the ΛCDM scenario, and Ω_m and w for the alternative scenario. We marginalize over the remaining cosmological parameters (h, Ω_b h², n_s, and τ for the first, and also σ₈ for the second model) by assigning uniform priors centered on their corresponding values in the PATCHY simulation cosmology (Planck Collaboration et al. 2014) and with a total amplitude equal to 10× the 68% uncertainty provided by the Planck2018 results (Planck Collaboration et al.2020). We tested the effect of using different priors, considering the constraints given by the Planck2013 cosmology (Planck Collaboration et al. 2014), finding no or nonsignificant variations (<1.5%) on both the predicted cosmological values and their relative precision.

Two fundamental quantities to compute at each step of the MCMC procedure are the Hubble parameter, H(z), and the matter power spectrum, P(k). We calculate the former as

$$\begin{eqnarray}&&H(z)={H}_{0}{\left[{{\rm{\Omega }}}_{{\rm{m}}}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{\mathrm{de}}{\left(1+z\right)}^{3(1+w)}\right]}^{1/2},\end{eqnarray} \tag{ 10 }$$

imposing w = −1 for the ΛCDM case, where Ω_de = Ω_Λ. The total matter power spectrum is instead evaluated with the approximations for the matter transfer function of Eisenstein & Hu (1999). All of the cosmological calculations (e.g., the 2PCF modeling reported in Appendix A), the Bayesian analyses, and the data catalog manipulations (including the void cleaning procedure described in Section 3.2) described in this work are performed in the numerical environment provided by the CosmoBolognaLib ¹³ (Marulli et al. 2016), a large set of free software C++/Python libraries composed of constantly growing and improving codes for cosmological analyses.

As discussed in Section 2, to model the size function of cosmic voids identified in a biased tracer field, we need to estimate the effective bias of our mass tracers. For this purpose, we measure the redshift-space anisotropic 2PCF of the analyzed galaxies by means of the popular Landy & Szalay (1993) estimator. We then model the corresponding multipole moments with the theoretical prescriptions presented in Taruya et al. (2010), which take into account both the nonlinear gravitational clustering, redshift-space distortions, and galaxy bias. We refer the reader to Appendix A for a detailed description of the 2PCF modeling performed in this work.

We apply the galaxy clustering modeling procedure to the 100 PATCHY mock galaxy catalogs obtaining the correlation matrix of the 2PCF multipoles and the values of b_eff σ₈ for each realization. By computing the median and the standard deviation of the 100 σ₈-dependent bias values of the PATCHY galaxy mocks, we attain b_eff σ₈ = 1.33 ± 0.02 and 1.26 ± 0.02, respectively for the two redshift bins considered in this analysis, i.e., 0.2 < z ≤ 0.45 and 0.45 < z < 0.65 (see Section 3.1). The corresponding values of the effective bias, achievable in this case knowing the true cosmology of the MultiDark PATCHY simulations, are b_eff = 1.94 ± 0.03 and 2.02 ± 0.03.

We use the values of b_eff σ₈ to calibrate the linear function ${ \mathcal F }$ reported in Equation (7) by means of a Bayesian analysis of the number counts of PATCHY cleaned voids. To this end, we measure the median number counts of the 100 PATCHY samples of voids processed by the cleaning algorithm, fixing the minimum void radius at 2.5× the mean galaxy separation (MGS) computed at the mean redshift of the galaxies in the low- and high-redshift subsamples considered. This choice is in agreement with previous works (Contarini et al. 2019, 2021, 2022; Ronconi et al. 2019) and leads to the selection of voids with approximately R_eff > 35 h⁻¹ Mpc. A detailed discussion of the impact of this choice is presented in Appendix C.

Our goal is to fit both the median PATCHY void count data sets (i.e., featuring the two selected redshift ranges) to calibrate the free parameters, C_slope and C_offset, of the extended Vdn model described in Section 2.2. As for B_slope and B_offset, these coefficients are meant to take into account the galaxy bias variation with redshift, b_eff(z), so our parameterization is expected to hold in the full range of redshift covered by the survey. Nevertheless we notice that, if the galaxy samples extracted at different redshifts were characterized by incompatible intrinsic properties (e.g., different flux or mass limits), this relation is in fact expected to vary (Contarini et al. 2019), eventually deviating from linearity. In our case, the galaxy sample selection is restricted to two redshift bins only, so the linear relation calibrated encompasses any potential variation in the intrinsic properties of the sample.

We remind the reader that the calibrated values, as well as the relative covariance relation for these parameters, are expected to be different from those found in previous papers (i.e., Contarini et al. 2019, 2021, 2022) for B_slope and B_offset. In these works indeed, ${ \mathcal F }$ was parameterized by means of b_eff only; therefore, the first coefficient of the linear function is positive, because of the increasing dependence of b_eff with redshift. ¹⁴ Conversely, the slope of the linear function ${ \mathcal F }$ is expected to be negative in our analysis given its dependence on b_eff σ₈, which is instead decreasing with redshift for the selected sample of galaxies.

We perform the Bayesian analysis starting by assuming the likelihood form expressed in Equation (9) and computing the covariance between our 100 PATCHY number count measures of the cleaned voids samples. We report the void counts covariance matrix, C _ij , normalized by its diagonal component, $\sqrt{{{\mathsf{C}}}_{{ii}}{{\mathsf{C}}}_{{jj}}}$, in the left panel of Figure 3. We separate in this plot the low- and high-redshift ranges considered in our analysis, characterized by six and five logarithmic radius bins, respectively. We notice that the covariance between the void radii at the same redshift and between the two redshift bins exhibits negligible off-diagonal terms (as suggested by Kreisch et al. 2022). We tested the effect of modeling the void counts with a Gaussian likelihood including the full covariance matrix presented here, and we found no statistically relevant change in the results, confirming the low impact of the off-diagonal terms in this analysis. We finally note that the supersample covariance effects are expected to be small on the void size function, so its covariance can be reliably approximated without considering the presence of supersurvey modes (Bayer et al. 2022).

Zoom In Zoom Out Reset image size

To model the average void number counts extracted from the PATCHY mock catalogs and to calibrate the void size function nuisance parameters, we applied the extended Vdn model to the low- and high-redshift bins analyzed in this work. We consider the median redshifts of the two galaxy samples, assigning to each the corresponding value of the σ₈-dependent bias as a Gaussian prior, with mean and standard deviation given by the 2PCF analysis. Then we set wide uniform priors for C_slope and C_offset, and sample the posterior distribution of these parameters, assuming the fiducial cosmology of the MultiDark PATCHY simulations.

We show the results of the nuisance parameter calibration in the right panel of Figure 3, where we report the 68% and 95% confidence contours on the parameters C_slope and C_offset. The values of the maximum posterior probability are reported at the top of the two panels representing the associated 1D posterior distribution projections. We notice that C_slope and C_offset are strongly anticorrelated, and that their posterior distribution is asymmetric. In particular, the skewness of the distribution is due to the physical requirement of ${\delta }_{{\rm{v}},\mathrm{DM}}^{\mathrm{NL}}\gt -1$ for the underdensity threshold, resulting from the Vdn model reparametrization.

In Figure 4 we report the calibrated Vdn model derived from the fitting procedure, together with the average void counts measured from the PATCHY mock catalogs. We show in the two panels of this figure the results for the low- and high-redshift bins analyzed, 0.2 < z ≤ 0.45 and 0.45 < z < 0.65. The markers indicate the number counts of cleaned voids, while the error bars are computed as the standard deviation of the 100 number count measures in each radius bin. We highlight that this uncertainty on the void counts is similar, although larger than the Poissonian error (i.e., the rms of the void counts), which represents an underestimation of the true error associated with the counts. The line indicates the theoretical prediction from the extended Vdn model, dependent on the value of b_eff σ₈, characterizing the galaxies belonging to the considered redshift subsample. We report with a shaded band the 68% uncertainty on the void size function model, associated with the constraints on the parameters C_slope and C_offset, derived from the calibration on the MultiDark PATCHY simulations.

Zoom In Zoom Out Reset image size

The agreement between the measured void number counts and the theoretical predictions of the calibrated Vdn model proves the effectiveness of our methodology, which we exploit in this work following two different approaches. The first is to assume the exact values, i.e., without any associated uncertainty, for the parameters C_slope and C_offset obtained from the calibration of the model on the PATCHY mocks. This case corresponds to the ideal assumption of using volume-unlimited cosmological simulations perfectly reproducing the survey data that one wants to study or, alternatively, to have a size function model including a fully theoretical description of both the effects of redshift-space distortions and mass tracer bias. This approach has been explored from a theoretical point of view in recent works (e.g., Verza et al. 2022), but a complete and calibration-independent model has not yet been provided. Although optimistic, this first approach shows the full potential of the methodology and is useful to understand the limitations related to the calibration procedure; the analysis performed with fixed values of C_slope and C_offset will be labeled hereafter as fixed calibration.

The second approach consists of including the uncertainties derived from the calibration procedure: the entire posterior distribution of the parameters C_slope and C_offset, sampled using the PATCHY void number counts and reported in the right panel of Figure 3, is assumed as prior for the Bayesian analysis of the BOSS void size function. The analysis conducted with this methodology will be denoted as relaxed calibration, and provides more conservative cosmological constraints, due to the limits of the calibration procedure performed with the MultiDark PATCHY simulations. We refer the reader to Appendix C for further insights about the impact of the calibration parameter priors on the cosmological constraints from the void size function.

To validate the proposed modeling approaches, we derive mock cosmological constraints from the PATCHY void size function, with the aim of recovering the true cosmology of the simulation and evaluating the constraining power of the measured void number counts. We report the results of this test in Appendix B, where we prove the effectiveness of the model in retrieving the MultiDark PATCHY simulation input cosmological parameters within 68% errors, for both the assumed cosmological models (i.e., ΛCDM and wCDM). We also refer the interested reader to Appendix B for a discussion about the dependency of the void size function constraining power on the main characteristics of the analyzed data: the number density and volume coverage of the galaxy sample. As these two quantities differ between the PATCHY and BOSS catalogs, the precision of the cosmological constraints obtained from the corresponding void samples is also slightly different.

Now that the void size function model has been calibrated, we are ready for the cosmological exploitation of the void number counts measured from the BOSS galaxy catalogs. To prepare the void sample, we follow the same pipeline applied to the PATCHY mocks. First, we clean the void catalogs extracted with VIDE, assuming the galaxy number density distribution reported in the left panel of Figure 1 and approximated with a spline function, which is represented as a dashed line. During the cleaning operation, we consider the same values for the density threshold, ${\delta }_{{\rm{v}},\mathrm{tr}}^{\mathrm{NL}}=-0.7$, and for the void size selection, R > 30 h⁻¹ Mpc. Then we extract from this sample the number of voids as a function of their radius for the two redshift ranges considered in this analysis (see Section 3.1). We finally compute the uncertainties on the binned void size function measures through a volume-based rescaling of the errors estimated for the PATCHY void number counts, as explained in Section 3.1 to describe the right panel of Figure 1.

These data sets, obtained with the fiducial Planck2013 cosmology (Planck Collaboration et al. 2014), need to be compared with the corresponding theoretical model of the void size function evaluated with the calibrated parameters C_slope and C_offset. To do this, we measure the value of the σ₈-dependent bias, considering the BOSS galaxies belonging to the two selected redshift subsamples, and the 2PCF covariance matrix built by means of PATCHY mocks (see Sections 3.4 and Appendix B): we find b_eff σ₈ = 1.36 ± 0.05 for the interval 0.2 < z ≤ 0.45 and b_eff σ₈ = 1.28 ± 0.06 for 0.45 < z < 0.65. Thereafter, we use these values to compute the theoretical void size function by means of the extended Vdn model calibrated in Section 3.4, and we rely on the effective volume covered by the void finder applied to the BOSS galaxies to convert the predicted void number density into number counts, directly comparable with our observed data sets.

At this point we perform a Bayesian analysis of the measured void number counts by assuming the likelihood function reported in Equation (9) and testing the two cosmological models outlined in Section 3.3. For both the considered cosmological scenarios, i.e., ΛCDM and wCDM, we follow the two approaches presented in Section 3.4, which depend on the constraints assumed from the model calibration with PATCHY mock catalogs, i.e., the fixed calibration and the relaxed calibration methodologies.

In Figure 5 we show a comparison between the measured size function of voids identified in the distribution of BOSS DR12 galaxies and the best-fit model computed with the MCMC sampling technique. We report the results for the ΛCDM cosmological model, obtained treating the parameters C_slope and C_offset of the extended Vdn model as free, but restrained by a prior distribution determined from the PATCHY calibration constraints (i.e., the relaxed calibration case). We recall that in this specific analysis, the free cosmological parameters of the model are Ω_m and σ₈, while the remaining ones are bound to Planck2018 constraints (see Section 3.3 for a more detailed description). The MCMC technique used to sample the posterior distribution of all of these parameters provides us with the best-fit value and the corresponding uncertainty for the void size function model, represented in Figure 5 as a solid line surrounded by a shaded band covering the 68% confidence region.

Zoom In Zoom Out Reset image size

We point out the good agreement between the measured void number counts and the extended Vdn model in Figure 5, for both the low- (left panel) and high- (right panel) redshift ranges analyzed, featuring three and four independent radius bins, respectively. We highlight that the void number count binning applied here is the same as that used for the PATCHY analysis. In this case, however, the void counts are not derived from an averaging procedure over different realizations and therefore cannot be lower than unity. This fact, combined with the statistical rarity of large voids, leads to more restricted void size ranges with respect to those presented in Figure 4. Moreover, analogously to what we found for the PATCHY mocks, voids are more abundant in the high-redshift subsample, given the greater volume covered by the galaxies in the light-cone.

The data-model agreement obtained for the remaining modeling approach (i.e., the fixed calibration case) and for the alternative cosmological scenario considered (i.e., wCDM, assuming both types of calibration parameter priors) is similar to the results presented, which is why we do not explicitly show them here. Nevertheless, the cosmological constraints corresponding to each case will be described in the following section.

Figure 6 shows the constraints resulting from both void size function calibration approaches presented in this work. In the left panel we show the results for the ΛCDM scenario, while in the right we show those for the alternative wCDM model. We represent in different colors the 68% and 95% confidence contours for the two cases labeled fixed and relaxed calibration. We note that when relaxing the calibration constraints, the corresponding confidence contours slightly enlarge as a consequence of the contribution given by the variation of the parameters C_slope and C_offset. The main effect of the relaxed calibration approach is to loosen the constraining power on the parameters σ₈ and w, for the ΛCDM and wCDM scenarios, respectively. On the other hand, as visible from the 1D marginalized posterior distributions reported in the subpanels of Figure 6, the corresponding impact on Ω_m is statistically negligible.

Zoom In Zoom Out Reset image size

One of the most important features to notice is the orientation of the Ω_m–σ₈ confidence contour reported in the left panel of Figure 6. First, this degeneracy direction is in agreement with the results of other works that analyzed the void size function in simulated catalogs (Bayer et al. 2021; Kreisch et al. 2022). Second, as demonstrated by Contarini et al. (2022) and Pelliciari et al. (2023), it is orthogonal to that of other standard cosmological probes, such as galaxy clustering, weak lensing, and cluster number counts. The latter is a fundamental characteristic to be exploited in order to maximize the constraining power given by probe combinations, thanks also to the negligible correlation of the void size function with standard statistics of overdense objects. We provide some insights on the physical origin of the orthogonality of these constraints in Appendix D.

We report the maximum posterior values and the relative 68% errors derived from the BOSS void size function analysis in Table 1. As we discuss in more detail in S. Contarini et al. (2023, in preparation), the constraints we derived are compatible with the Planck2018 results (Planck Collaboration et al. 2020). Assuming the flat ΛCDM framework, we obtain a value of Ω_m = 0.285 (Ω_m = 0.293) for the relaxed (fixed) calibration approach, with a corresponding precision of about 21%. In the same model, we have σ₈ = 0.793 (σ₈ = 0.839) for the relaxed (fixed) calibration approach, with a precision of 10% (7%).

Table 1. Maximum Posterior Values and 68% Uncertainties Associated with the Confidence Contours Reported in Figure 6

Model Calibration	ΛCDM				wCDM
	Ωm	σ8	Cslope	Coffset	Ωm	w	Cslope	Coffset
Relaxed	${0.285}_{-0.057}^{+0.063}$	${0.793}_{-0.076}^{+0.090}$	$-{2.02}_{-0.61}^{+0.66}$	${4.74}_{-0.82}^{+0.87}$	${0.293}_{-0.047}^{+0.055}$	−1.11 ± 0.24	$-{2.07}_{-0.73}^{+0.70}$	${4.90}_{-0.97}^{+1.00}$
Fixed	${0.293}_{-0.059}^{+0.067}$	${0.839}_{-0.058}^{+0.055}$	−1.62	4.30	${0.299}_{-0.046}^{+0.051}$	$-{0.97}_{-0.16}^{+0.15}$	−1.62	4.30

Notes. We present the results for the two cosmological scenarios analyzed in this work, i.e., ΛCDM and wCDM, in two separate sets of columns, grouping together the freely varied parameters of the corresponding model. The first row of the table presents the results obtained by modeling the size function of BOSS DR12 voids using the relaxed calibration method, while the second row reports the constraints when adopting the fixed calibration approach (see Section 3.4 for further details).

Download table as:  ASCIITypeset image

As already mentioned, we expect a significant improvement of the precision on these parameters through the joint analysis of void counts with different cosmological probes. This is investigated in a companion paper (S. Contarini et al. 2023, in preparation), where we present a more specific and detailed analysis of the contribution given by cosmic void counts on two of the main current issues in cosmology: the tension between late and early probes on the Hubble constant H₀,and the matter clustering strength parameterized by ${S}_{8}\equiv {\sigma }_{8}\sqrt{{{\rm{\Omega }}}_{{\rm{m}}}/0.3}$ (see, e.g., Di Valentino et al. 2021a,2021b, for a review).

For what concerns the flat wCDM scenario, we obtain Ω_m = 0.293 (Ω_m = 0.299), with a corresponding precision of 17% (16%) for the relaxed (fixed) calibration strategy. Then we compute w = −1.11 (w = − 0.97) for the relaxed (fixed) calibration approach, with a precision of 22% (16%). For both cases, the results are compatible with the constant w = −1, so in agreement with the benchmark ΛCDM model.

We perform here the analysis on the PATCHY mocks by increasing the statistical errors of the void number counts. Specifically, during the calibration procedure, we associate to void counts Poissonian errors multiplied by factors of 2, 3, and 4. This approach is aimed at reproducing the effect of calibrating the void size function using mock catalogs with volumes 4×, 9×, and 16× smaller, respectively.

In the upper panel of Figure C1 we show the 2D posterior distribution of the parameters C_slope and C_offset, for the three analyzed cases. We note that the enlargement of the confidence contours is mainly restricted to the degeneracy direction of the two calibration parameters. In the lower panels of Figure C1 we show instead the test constraints obtained for the parameter planes Ω_m–σ₈ (ΛCDM scenario) and Ω_m–w (wCDM scenario). In all of the presented cases, the true values of the simulation cosmological parameters remain in the 68% confidence region of the posterior distributions. As expected from the comparison between the relaxed calibration and fixed calibration methods (see Figure B1), the impact of the nuisance parameter prior on Ω_m is almost negligible, especially in the ΛCDM case. In fact, the relative error on Ω_m is stable around 14%. The effect on the constraints on σ₈ and w is instead more important because of their degeneracy with the calibration parameters. The relative error on σ₈ increases from 9%–23%, when increasing the calibration errors from 2–4× the Poissonian error. Analogously, the relative error on w increases from 18%–25%. However, we underline that the effects shown in this test are derived from extreme calibration scenarios.

Zoom In Zoom Out Reset image size

In this section we repeat the whole analysis, assuming a fiducial cosmology that is different from that of Planck2013. We run the void finder and the cleaning procedure on the 100 PATCHY catalogs, considering a value of Ω_m either decreased or increased by 3× the standard deviation reported in Planck Collaboration et al. (2020). Specifically, for this test we assume Ω_m = 0.279 (−3σ) and Ω_m = 0.335 (+3σ), modifying the value of Ω_de to maintain Universe flatness. Moreover, we highlight that the accepted void radius range is recomputed during the analysis according to the value of MGS(z) in the considered cosmology.

We show in Figure C2 the 2D posterior distributions obtained for the parameter spaces Ω_m–σ₈ and Ω_m–w. The two confidence contours are statistically compatible with each other, and the true cosmological values of the PATCHY simulations are recovered within the 68% uncertainty in both cases.

Zoom In Zoom Out Reset image size

We report in Table C1 the corresponding marginalized constraints computed, for the ΛCDM and wCDM scenarios. We note that, especially when exploring fiducial Ω_m values lower than the true one (see the Ω_m = 0.279 case), the underestimation of these cosmological parameters starts to become statistically significant. This issue will be analyzed in detail in future works.

Table C1. Maximum Posterior Values and 68% Uncertainties of the Cosmological Constraints Computed Assuming Ω_m = 0.279 and Ω_m = 0.335 as Fiducial

ΛCDM	Ωm	σ8
Fiducial Ωm = 0.279	${0.265}_{-0.034}^{+0.041}$	${0.775}_{-0.064}^{+0.089}$
Fiducial Ωm = 0.335	${0.270}_{-0.034}^{+0.042}$	${0.786}_{-0.067}^{+0.084}$
wCDM	Ωm	w
Fiducial Ωm = 0.279	${0.278}_{-0.033}^{+0.035}$	$-{1.18}_{-0.21}^{+0.22}$
Fiducial Ωm = 0.335	${0.284}_{-0.036}^{+0.035}$	$-{1.17}_{-0.20}^{+0.24}$

Note. We report the results for both the ΛCDM (upper rows) and wCDM (lower rows) scenarios. The associated confidence contours are shown in Figure C2.

Download table as:  ASCIITypeset image

As a final test on systematic errors, we focus here on the void radius selection. The choice regarding the minimum acceptable void radius holds significant importance, as it directly impacts the strength of the resulting constraints: the greater the number of small voids included in the analysis, the more stringent the constraints derived from the void number counts. However, it is equally crucial to exclude the spatial scales that are influenced by void count incompleteness. In fact, the resolution of the tracer catalog affects the identification of small voids, which may be lost because of the Poisson noise.

In this work, we restricted our cleaned void sample to radii larger than 2.5× the MGS of the galaxy catalog (see Sections 3.1 and 3.4). This leads to a minimum void radius of about ${R}_{\min }=36.37$ and 35.13 h⁻¹ Mpc, for the redshift bins 0.2 < z ≤ 0.45 and 0.45 < z < 0.65, respectively. For simplicity, we apply the same radius cut to both redshift bins, without relying on the MGS value. Here, we consider three different radii selections: one slightly lower than that used in the main analysis of this work, i.e., ${R}_{\min }=35\ {h}^{-1}\ \mathrm{Mpc}$, and two higher, i.e., ${R}_{\min }=37.5$ and 40 h⁻¹ Mpc.

In Figure C3 we show the 2D posterior distributions computed for the ΛCDM and the wCDM scenarios by applying the presented radius selections. As expected, the cosmological constraints become weaker as the void minimum radius becomes larger. In the ΛCDM scenario, moving from the lowest to the highest radius cut, the relative error increases by factors of 1.9 and 1.2, for Ω_m and σ₈, respectively. When considering the wCDM scenario, the increase is instead by a factor 1.6, for both the parameters Ω_m and w. However, it is important to emphasize that, in all three cases, the true cosmological parameters of the PATCHY simulations are within the 68% confidence region.

Zoom In Zoom Out Reset image size