© The Authors 2025

1. Introduction

For over 40 years, it has been observed that the brightest galaxies are generally found in dense regions and that most of cosmic space is devoid of these types of galaxies. This distribution is amply understood as a natural evolution of density fluctuations in matter that have been progressively amplified by gravitational instability (Giovanelli 2010). In this framework, it has been demonstrated that initially under-dense regions grow in size, while highly dense areas end up collapsing under their own gravity (Sheth & Van De Weygaert 2004).

Over the last century, there has been a significant focus on studying the over-dense regions of the Universe (e.g. Kiang & Saslaw 1969; Bahcall 1977; Kaiser 1987; Einasto et al. 1994; Holder et al. 2001; Yang et al. 2005; Li et al. 2006; Gao & White 2007; Zentner et al. 2019; Dong-Páez et al. 2024), while under-dense regions only started attracting interest more recently (e.g. Hoffman & Shaham 1982; Zeldovich et al. 1982; Blumenthal et al. 1992; Little & Weinberg 1994; Goldwirth et al. 1995; Sheth & Van De Weygaert 2004; Li et al. 2012; Pisani et al. 2014; Pollina et al. 2016; Achitouv 2019; Chan et al. 2019a; Rodríguez-Medrano et al. 2024; Curtis et al. 2024). Studying these regions (voids) is very useful as they possess unique characteristics that make them important probes for cosmological studies and the physics of galaxy formation. They are useful, for example, in:

• Constraining the equation of state for dark energy (e.g. Lee & Park 2009; Biswas et al. 2010; Sutter et al. 2014; Contarini et al. 2022),

• Studying modified gravity (e.g. Martino & Sheth 2009; Clampitt et al. 2013; Voivodic et al. 2017; Falck et al. 2017; Perico et al. 2019; Contarini et al. 2021; Mauland et al. 2023)

• Constraining cosmological models (e.g. Ryden 1995; Benson et al. 2003; Croton et al. 2005; Lavaux & Wandelt 2010),

• Constraining cosmological parameters based on their statistics (e.g. Betancort-Rijo et al. 2009; Nadathur 2016; Hamaus et al. 2020; Aubert et al. 2022; Contarini et al. 2022, 2023),

• Testing the primordial non-Gaussianities (e.g. Song & Lee 2009; Chongchitnan & Silk 2010; Chan et al. 2019b),

• Studying massive neutrinos and dark energy (Verza et al. 2019; Contarini et al. 2022; Verza et al. 2023; Thiele et al. 2024).

In fact, any cosmological parameter could be constrained, in principle, from void statistics. Specifically, the abundance of voids with radii larger than r, N_v(r) depends on the normalisation of the linear spectrum of density fluctuations (σ₈). The shape of this function for large values of r depends on Γ = Ω_ch, where h = H₀/100 km s⁻¹ Mpc⁻¹ (see Betancort-Rijo et al. 2009 for details).

Constraining σ₈ and Γ = Ω_ch is especially interesting as there some statistically significant tension does exist for the cosmological parameter S${}_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ between the Planck experiment (Planck Collaboration VI 2020), which measures the cosmic microwave background (CMB) anisotropies (S₈ = 0.834 ± 0.0161), and other low-redshift cosmological probes, such as weak gravitational lensing, where a value of S${}_8=0.{776}_{-0.030}^{+0.032}$ was obtained by the Kilo Degree Survey (KiDS-1000) (Li et al. 2023). This tension is known as the S₈ tension. (see Di Valentino et al. 2021 for more details).

Early studies of cosmic voids have been limited by the tiny areas inspected by the available surveys. Nevertheless, with the emergence of more expansive redshift surveys, such as Two-Degree Field Galaxy Redshift Survey (2dFGRS) (Colless et al. 2001, 2003), Sloan Digital Sky Survey (SDSS) (York et al. 2000) and the SDSS Baryon Oscillation Spectroscopic Survey (BOSS) (Dawson et al. 2013), alongside improved resolution in cosmological simulations and enhanced analytical methodologies, we now have the capability to derive precise statistical insights concerning voids (e.g. Tikhonov 2006, 2007; Ceccarelli et al. 2006; Von Benda-Beckmann & Müller 2008; Patiri et al. 2012; Hamaus et al. 2020; Douglass et al. 2023; Contarini et al. 2023).

However, despite the longstanding presence of the concept of cosmic voids, there is no universally accepted definition for what constitutes a void. The term “void” can encompass a range of disparate entities, depending on the data employed and the objectives of the analysis. For example, voids can be defined as under-dense regions based on the smoothed dark matter (or halo/galaxy) density field (Colberg et al. 2005), as gravitationally expanding regions based on the dynamics of the dark matter density field (Hoffman et al. 2012; Elyiv et al. 2015) or as empty spatial regions among discrete tracers (El-Ad & Piran 1997; Aikio & Mähönen 1998; Hoyle & Vogeley 2002).

The choice of a simple definition of voids, in particular, their definition as empty spheres, is convenient for statistical studies of galaxy voids. In this paper, we define voids as maximal non-overlapping spheres empty of objects with mass (or luminosity) above a given value. With this definition, it is clear that voids are not empty, as there can be low luminosity galaxies (or low mass halos) inside them.

The aim of this paper is to infer the values of σ₈, H₀, and Ω_m from SDSS redshift survey using the theoretical framework of voids statistics developed in Betancort-Rijo (1990). This has been already done in previous articles in different ways. For example, in Sahlén et al. (2016) galaxy cluster and void abundances are combined using extreme-value statistics on a large cluster and a void. In this way, they have been able to obtain a value of σ₈ = 0.95 ± 0.21 for a flat ΛCDM universe. In Hamaus et al. (2016), the authors constrained Ω_m = 0.281 ± 0.031 by studying void dynamics in SDSS. In Woodfinden et al. (2023) they also used SDSS survey to constrain ${\mathrm{\Omega }}_{\mathrm{m}}=0.{391}_{-0.021}^{+0.028}$ by measuring the void-galaxy and galaxy-galaxy clustering. In Contarini et al. (2023) they modelled the void size function (Press & Schechter 1974; Sheth & Van De Weygaert 2004) by means of an extension of the popular volume-conserving model (Jennings et al. 2013), which is based on two additional nuisance parameters. Applying a Bayesian analysis, they obtained a value of ${\sigma }_8=0.{79}_{-0.08}^{+0.09}$ for the BOSS DR12 survey; in a posterior work (Contarini et al. 2024), they constrained S${}_8=0.{813}_{-0.068}^{+0.093}$ and H${}_0=67.3_{-9.1}^{+10}$ for the same redshift survey. However, in these works, the uncertainties obtained for σ₈ and H₀ are very large and the constrained values are compatible with high- and low-redshift cosmological probes. In addition, none of these studies have sought to combine their results with the CMB or with weak lensing results. Therefore, in this work, we propose using a original theoretical framework, presented for the first time in Betancort-Rijo et al. (2009), which we have calibrated and tested on different halo simulation boxes in this work. We inferred σ₈, Ω_m, and H₀ from SDSS galaxies with redshifts of 0.02 < z < 0.132 and absolute magnitudes in the r-band of M_r < −20.5. We combined these results with KiDS-1000 and the Dark Energy Survey (DESY3) weak lensing results.

This paper is structured as follows. In Section 2, we describe the redshift survey as well as the simulations used in this work. In Section 3, we give a detailed explanation of how our ‘Void Finder’ works and give the relevant statistics of voids obtained for the observational catalogue, comparing it with the result of Uchuu-SDSS lightcones. In Section 4, we introduce the most important concepts and equations used in the theoretical framework used to calculate the abundance of voids larger than r and the void probability function (VPF). We demonstrate that the theoretical framework successfully predicts these two void statistics for halo simulation boxes with different σ₈ values in real and redshift space, and we show the results for SDSS and the Uchuu-SDSS lightcones. In Section 5, we explain the statistical test we used to infer σ₈, Ω_m, and h. In Section 6, we show the constrained values for the sample of Uchuu-SDSS lightcones and SDSS redshift survey we used in this work and we combined our results with KiDS-1000 results (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023). In Section 7, we compare our constrained values of σ₈, Ω_m, and H₀ with those obtained in other works where voids are also used. Finally, in Section 8, we summarise the most important results obtained in this work.

2. Data and mocks

The aim of this work is to infer σ₈, Ω_m, and H₀ values from Sloan Digital Sky Survey. However, to make sure that the theoretical equations correctly reproduce the void functions of this redshift survey, we also used four halo simulation boxes with different σ₈ values (see Table 1), one simulation galaxy box and 32 lightcones. In this section, we introduce the redshift survey and mocks.

2.1. SDSS

We used the seventh release (Abazajian et al. 2009) of Sloan Digital Sky Survey (SDSS DR7) (York et al. 2000), which includes 11 663 deg² of CCD imaging data in five photometric bands for 357 million distant objects. The catalogue also completed spectroscopy over 9380 deg². In total, there are 1.6 million spectra, including 930 000 galaxies, 120 000 quasars, and 460 000 stars.

However, in this work we use a subcatalogue of SDSS: only galaxies from the northern regions with completeness greater than 90% are selected. Therefore, the effective area of this subcatalogue is 6511 deg² and contains around 497 000 galaxies with redshifts between z ∈ (0, 0.5). In this sample, ∼6% of targeted galaxies lack a spectroscopically measured redshift due to fibre collisions. Thus, a nearest neighbour correction was applied to these galaxies, assigning to them the redshift of the galaxy they have collided with (Dong-Páez et al. 2024).

Additionally, we imposed further cuts on absolute magnitude and redshift. We constructed a volume-limited sample by keeping only galaxies brighter than the Milky Way-like galaxies (M_r < −20.5, where M_r is the absolute magnitude in r-band) and with z ≤ 0.132. We additionally imposed z ≥ 0.02, so that we did not take into account any nearby galaxies highly affected by peculiar velocities. The physical volume of this sample is V = 41.67 × 10⁶ h⁻³ Mpc³. There is a total of 112 496 galaxies that fulfill these restrictions, so the number density of the galaxies is n_g = 2.838 × 10⁻³ h³ Mpc⁻³.

2.2. Mocks

We used the Uchuu simulation, which was produced using the TreePM code GreeM¹ (Ishiyama et al. 2009, 2012) on the supercomputer ATERUI II at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. The number of dark matter particles was 12800³ with a mass resolution of 3.27 × 10⁸ h⁻¹ M_⊙ in a box with a side length of 2000 h⁻¹ Mpc. A total of 50 halo catalogues (snapshots) were created, covering the redshift range from 0 to 14 (Ishiyama et al. 2021). The halos were subsequently identified using the RockStar halo/subhalo finder² (Behroozi et al. 2013a) and merger trees were generated using the consistent trees code³ (Behroozi et al. 2013b). These simulations adopted the cosmological parameters from Planck 2015 (see Table 1) (Planck Collaboration XIII 2016). All this data is publicly available and accessible in the Skies & Universes database⁴, including galaxy catalogues constructed using various methods (Aung et al. 2023; Oogi et al. 2023; Gkogkou et al. 2023; Prada et al. 2023; Ereza et al. 2024; Dong-Páez et al. 2024).

Apart from Uchuu, we use three more simulations which are being used for the first time in this work. One of these simulations have Planck 2018 cosmological parameters (Planck Collaboration VI 2020) P18 and the other two have Planck 2018 parameters, but with σ₈ = 0.75 (P18LowS8) σ₈ = 0.65 (P18VeryLowS8). Additionally, P18 and P18LowS8 have exactly the same simulation properties, except for the value of σ₈. P18VeryLowS8 has different value of σ₈ and different mass resolution. All these details can be seen in Table 1. These three simulations were produced using the TreePM code GreeM on the supercomputer Fugaku at the RIKEN Center for Computational Science. As well as the Uchuu simulation, we generated initial conditions for these three simulations with 2nd order Lagrangian perturbation theory (Crocce et al. 2006) by 2LPTIC code⁵. The initial conditions are identical across the three simulations, enabling us to compare them without cosmic variance. The total of 70 halo catalogues were created, where the redshift list is the same with the Uchuu at z < 6. N-body data, including the halo catalogues and the merger trees are also available in the Skies & Universes database.

The halo mass function (HMF) of each simulation can be seen in Figure 1, where M_vir, all is the halo mass within virial radius including unbound particles. In this figure, it can be seen that VeryLow has lower mass resolution, as the HMF decreases considerably for small masses. It is also shown that (especially for large masses) the HMF power decreases proportionally with σ₈. Additionally, it can be seen that there is no significant difference between P18, Low and VeryLow halo mass functions for low virial masses (M_vir, all < 10¹² h⁻¹ M_⊙), but there is a significant difference for large virial masses (M_vir, all > 10¹² h⁻¹ M_⊙).

Fig. 1. Halo mass function multiplied by the mean halo mass within virial radius (including unbound particles) of each bin of the four box catalogues used in this work. Shaded regions represent poissonian errors.

We also used a galaxy box with SDSS properties and 32 Uchuu-SDSS lightcones constructed with these galaxy boxes at different snapshots. These products are presented and extensively detailed in Dong-Páez et al. (2024). The galaxy box has been generated from the Uchuu simulation using the subhalo abundance matching (SHAM) method. Lightcones with SDSS properties have been constructed from the Uchuu-SDSS galaxy boxes, using a total of six snapshots between z = 0 and z = 0.5, which are separated in redshift by approximately 0.1 (z_snap = 0, 0.093, 0.19, 0.3, 0.43, and 0.49). We refer to Dong-Páez et al. (2024) for more information about the methodology followed to construct the lightcones.

For the purpose of studying void statistics in simulation boxes, we select those snapshots corresponding to a redshift of z ∼ 0.092 (which is the snapshot with the closest redshift to SDSS median redshift). Moreover, we did not use all the objects (halos or galaxies) from this boxes, but we select a number density of $\overline{n}=3\times {10}^{-3}{h}^3$ Mpc⁻³ for each box. This can be reached removing halos less massive than M_vir, all = 1.626×10¹² M_⊙/h (Uchuu), 1.642 × 10¹² M_⊙/h (P18), 1.595 × 10¹² M_⊙/h (P18LowS8), and 1.457 × 10¹² M_⊙/h (P18VeryLow), along with galaxies that are less bright than M_r < −20.5 for Uchuu galaxy box.

The number density of the boxes, $\overline{n}=3\times {10}^{-3}{h}^3$ Mpc⁻³, does not exactly match the number density of SDSS or Uchuu-SDSS due to the incompleteness of the latter. However, by applying the same cut in M_r, we could ensure that the objects in the boxes are consistent with those in Uchuu-SDSS, except for those “lost” due to incompleteness. Additionally, halo boxes were used to calibrate the theoretical framework and test its ability to successfully predict void functions for different values of σ₈.

3. Statistics of voids in the SDSS

The next important step we needed to take once the sample of halos or galaxies we are going to use is well defined is to then identify voids, defined as maximum non-overlapping spheres. There are a number of previous works that have used the same definition of voids and have developed their own void finder, such as (Jennings et al. 2013; Ronconi & Marulli 2017; Contarini et al. 2024; Zhao et al. 2016). In this work, we use the same procedure as in Zhao et al. (2016), namely:

1. We performed a Delaunay triangulation (with periodic conditions for simulation boxes),

2. Delaunay triangulation provides us with a list of spheres that fulfill the Delaunay condition: the Delaunay triangulation of a set of points, p_i, ensures that no point, p_i, lies within the circumcircle of any triangle in the triangulation. However, these spheres can overlap, so we write an additional code that find voids among these spheres.

In order to calculate the Delaunay triangulation, we use the public code called Delaunay trIangulation Void findEr (DIVE⁶), which is introduced in Zhao et al. (2016). This code uses The Computational Geometry Algorithms Library⁷, and, as mentioned above, the output of this code is a file that contains the positions in space of the centres of the spheres that satisfy the Delaunay condition, as well as their radii. However, these spheres are not voids, but candidates to be voids.

To find voids (i.e. maximal non-overlapping spheres) among these spheres, an additional code needs to be developed. This code must check if two spheres overlap and, in case they do, keep the largest one as a void (see Appendix A for a more detailed explanation of the algorithm).

In Figure 2 the voids larger than r > 9 h⁻¹ Mpc found in the P18 halo simulation box with number density 3 × 10⁻³ h³ Mpc⁻³ are shown in a region of the box. It can be seen that each void is defined by four galaxies laying in its surface (some voids in Figure 2 have fewer than four galaxies due to being outside the plotted region) and we could ensure that voids do not overlap.

Fig. 2. Voids with r > 9 h−1 Mpc (spheres) found in a region of P18 box catalogue (black and red points) with number density 3 × 10−3 h3 Mpc−3. Points that define the voids (i.e. those lying in their surface) are highlighted with a red circle. The volume of the box is 803 h−3 Mpc3.

However, for the galaxy lightcones, one more step is required. This step involves considering the incompleteness in stellar mass, which causes many spurious voids to appear. These spurious voids correspond to regions in the catalogue with very low completeness, where galaxies could not be detected properly. Therefore, these false voids must be eliminated using an angular mask, such as the Healpix map (Gorski et al. 1999; Blanton et al. 2005; Swanson et al. 2008), characteristic of the redshift survey.

The procedure for removing these false voids is as follows: first, we generated points uniformly distributed within the volume of the void. Next, we projected these points in the angular plane, calculated the completeness of each point, and averaged the completeness of all points within the void to obtain an approximation of the completeness of the void. If this completeness is equal to or greater than 0.9, we labelled that void as a true void. Otherwise, we removed the false void.

Additionally, voids whose centres or part of their volume are outside the lightcone in the radial direction must also be removed, as the procedure explained above considers only the angular plane; however, voids may extend beyond the sample in the radial direction, so this extra check is necessary.

In Figure 3, the number density of voids larger than r, n_v(r), obtained for voids found in the sample of SDSS and Uchuu-SDSS (the mean of the 32 lightcones) considered in this work (galaxies with M_r < −20.5 and 0.02 < z < 0.132) are shown (in Table 2 the values are multiplied by the volume). It can be checked that Uchuu-SDSS statistics are compatible with observations for all radius bins, and both are compatible with the values predicted by the theory.

Fig. 3. Number density of voids larger than r, nv(r)c, obtained for voids found in the distribution of SDSS (black points) and Uchuu-SDSS galaxies (pink solid line) with Mr < −20.5 and 0.02 < z < 0.132 and predicted by the theoretical framework described in Section 4 with Planck 2015 parameters (black dashed line). Shaded pink region delimits the standard deviation (1σ) of the 32 Uchuu-SDSS lightcones.

The volume used in order to calculate the number density of voids larger than r in Figure 3 is:

$$\begin{array}{c}\hfill V(r)=\frac{N_v(r)(\mathrm{Uchuu}-\mathrm{SDSS})}{n_v(r)(\mathrm{Uchuu})},\end{array}$$(1)

where N_v(r)(Uchuu − SDSS) is the number of voids larger than r found in Uchuu-SDSS lightcones and n_v(r)(Uchuu) is the number density of voids larger than r found in Uchuu galaxy box. Although this quantity has units of volume, it is not a real volume. However, this quantity allow us to convert lightcone statistics into number densities. This volume takes into account all the effects that can impact the volume that a real void can occupy, such as the completeness (amongst others).

An important remark about Figure 3 is that if we want to compare the observed number density of voids larger than r obtained from SDSS (or Uchuu-SDSS) with that given by the theoretical framework presented in this work, we must take into account that the former suffers from incompleteness, as well as other effects such as border effects, while the latter does not. Therefore, we have to transform SDSS (and Uchuu-SDSS) void statistics as if it did not suffer from these effects. One way of doing this is using Uchuu galaxy box. Then, SDSS and Uchuu-SDSS number density of voids larger than r must be transformed as:

$$\begin{array}{c}\hfill n_v(r)\to n_v(r)(\mathrm{Uchuu})\frac{N_v(r)(\mathrm{SDSS})}{N_v(r)(\mathrm{Uchuu}-\mathrm{SDSS})},\end{array}$$(2)

where N_v(r)(Uchuu − SDSS) is the number of voids larger than r found in Uchuu-SDSS lightcones, while N_v(r)(Uchuu) is the same but found in Uchuu galaxy box.

4. Void statistics theoretical framework

The void size function of voids was predicted for the first time Sheth & Van De Weygaert (2004), where the authors used the same excursion-set approach used for the mass function of dark matter halos Press & Schechter (1974). This model is also based on the definition of voids as under-dense, non-overlapping spheres, which have gone through shell crossing (see Press & Schechter 1974 for more details). However, this model predicts that the fraction of the volume occupied by voids can be larger than the total volume of the universe. This was fixed by Jennings et al. (2013), where they introduce the volume-conserving model (Vdn model), where the total volume occupied by cosmic voids is conserved in the transition to the non-linear regime.

As demonstrated in Ronconi et al. (2019), the Vdn model is able to accurately predict the measured void size function of unbiased tracers, provided that the void catalogue is approximately cleaned from spurious voids and that the void radii are re-scaled to a fixed density threshold. This means that the developed void finder has to find voids with these requirements. However, in many of the works where this model is used (Contarini et al. 2019, 2022, 2023, 2024) a void finder based on Voronoi tesselations has been employed, which does not require voids to be spherical; then, a cleaning algorithm was applied to remove small voids and transform the remaining voids, so that they are non-overlapping maximal spheres, with a specific density contrast threshold.

Additionally, the Vdn model predicts the number density of voids found in the distribution of dark matter particle field. In Contarini et al. (2019), an extension of this model is proposed by introducing a parametrisation of the Vdn model’s non-linear under-density threshold as a function of the tracer effective bias, which has two free parameters that must be calibrated using simulations.

In this work, we propose a new model that only requires that voids be defined as maximal non-overlapping spheres. Thus, when using our void finder, we do not need to carry out a cleaning procedure. Additionally, the bias is completely determined by the model.

In this section, we present the theoretical framework used in this work, which was present for the first time in Betancort-Rijo et al. (2009), to calculate the number density of voids larger than r. First, we show the equations in real space and then we explain how the equations change in redshift space. We compare the results provided by the theoretical framework with the results obtained with the four halo simulations to check its validity for different values of σ₈ and in real and redshift space. The last step is to compare the model with SDSS void statistics.

4.1. Real space

The main void statistics we will study in this work is the number density of voids larger than r, ${\overline{n}}_v(r)$⁸, which can be related to void probability function, P₀(r) (White 1979) using the (re-calibrated) expression from Betancort-Rijo et al. (2009):

$$\begin{array}{c}\hfill {\overline{n}}_v(r)=\frac{\mu \mathcal{K}(r)}{V}{e}^{-\alpha \mathcal{K}(r)},\end{array}$$(3)

where μ = 0.588 and α = 1.671. These values have been calibrated using Uchuu simulation (see Appendix E) and

$$\begin{array}{c}\hfill \mathcal{K}(r)={[-\frac{1}{3}\frac{dln{P}_0(r)}{\mathit{dlnr}}]}^3{P}_0(r),\end{array}$$(4)

with

$$\begin{array}{c}\hfill V=\frac{4}{3}\pi r^3.\end{array}$$(5)

In Eq. (3), it is assumed that ${\overline{n}}_v(r)$ is an universal functional of P₀(r), so that independently of the clustering process underlying P₀(r), ${\overline{n}}_v(r)$ is related to it by an unique expression. However, it has been shown that for white noise, μ = 0.68. Thus, it is clear that this coefficient has a dependence on the clustering properties of the objects considered. However, we have found that for the simulations that we have used the quoted values of μ and α are valid. Therefore we shall use these values in all our considerations.

It is important to remark that Eq. (3) is only valid for 𝒦(r)≤0.46. For 𝒦(r) > 0.46, ${\overline{n}}_v(r)=0.313/V$. 𝒦(r), measures the rareness of the voids. The expression of P₀(r), which is the void probability function (VPF) is derived in Appendix B and in Betancort-Rijo et al. (2009), where it is shown that this expression arises from first principles. Additionally, in Appendix D, we show that the VPF successfully reproduces the VPF of the four simulation boxes used in this work.

Equations (B.11) and (B.12) were fit considering a cosmology with σ₈ = 0.9, so these expressions are only valid for values of σ₈ around this value (Rubiño-Martín et al. 2008). Moreover, it was demonstrated that small voids do not carry any cosmological information, so in this work, we only consider voids with r > 10 h⁻¹ Mpc. From Eqs. (3) and (4) we can see that if we know the VPF, we automatically know ${\overline{n}}_v(r)$; thus, we chose to first study the VPF and then calculate ${\overline{n}}_v(r)$ using Eq. (3).

We can define P₀(r) from a more general statistic, that is, P_n(r): the probability that a sphere of radius, r, placed at random within the distribution, contains n objects. If we assume a Poisson process, we can then write P_n(r) as (Layzer 1954):

$$\begin{array}{c}\hfill P_n(r)={\displaystyle {\int }_0^{\infty }}P(u)\frac{u^n}{n!}{e}^{-u}du,\end{array}$$(6)

where P(u) is the probability distribution for the integral of the probability density, u, within a randomly placed sphere. It is important to note that the dependence of P_n(r) with r is through u(r) (see Eqs. (B.6) and (B.7)).

An explicit computation of each term in P_n(r) can be consulted in Appendix B. There, it can be seen that P_n(r) depends on σ₈ (Eq. (B.16)), Γ = Ω_ch (Eqs. (B.11), (B.18), and (B.19)), and Ω_m (Eqs. (B.20) and (8)). However, the dependence on Ω_m is small.

We can see that so far our theoretical framework depends on four cosmological parameters: σ₈, Ω_c, H₀, and Ω_m. We note that the latter does not depend on Ω_Λ, as we impose a flat ΛCDM model. Also, the dependence of the theoretical framework with Ω_c and H₀ is only through Γ = Ω_ch. Additionally, the dependence on Ω_m is very weak, as it only enters through D(a). However, we can decrease the number of free parameters by relating Ω_m and Ω_b = Ω_m − Ω_c, where Ω_b is the baryonic matter density parameter. From Planck 2018 (Planck Collaboration VI 2020) we get α ≡ Ω_b/Ω_m = 0.157. Therefore, Ω_c = (1 − α)Ω_m and Ω_b = αΩ_m, so our final parameters are σ₈, Ω_m and H₀.

Before using the theoretical framework to constrain σ₈, Ω_m, and H₀ in the SDSS survey, it is important to check if the theoretical equations are accurate and have enough precision to recover the real values of σ₈, Ω_m, and H₀ from simulations with different values of the parameters. To do this, we used the halo simulation boxes that we have presented in Section 2.2, which have different values of σ₈.

The upper panel of Figure 4 gives the number density of voids larger than r predicted by the theoretical framework (continuous lines) and obtained by simulations (points) for the four halo simulation boxes in real space can be seen. In the bottom panel of the same figure, the ratio between simulations and theoretical framework is shown. Again, the agreement between the theoretical framework and the simulation values is good, especially for r between 12 and 18 h⁻¹ Mpc.

Fig. 4. Number density of voids larger than r in real space is shown for the theoretical framework (lines) and simulations (dots), shown in the left panel. The ratio between simulations and theoretical framework is shown in the bottom panels for the Uchuu, P18, Low, and VeryLow box catalogues with a number density of $\overline{n}=3\times {10}^{-3}{h}^3$ Mpc−3. The shaded region in the bottom panel delimits the region between 0.9 and 1.10 for the ratio.

4.2. Redshift space

The theoretical framework developed above is only valid for voids in real space. However, if we want to constrain σ₈ using surveys such as SDSS, this surveys provide galaxy positions in redshift space. In this space, the peculiar velocity of the galaxies is added to the velocity expansion of the Universe. This generates some distortions which result in elongated structures known as Fingers of God (see Hamilton 1998 for more details). If we want to calculate ${\overline{n}}_v(r)$ and P_n(r) in redshift space instead of real space, all we have to do is adjust r via:

$$\begin{array}{c}\hfill r=r^{\ast }\times {\{1+g\mathrm{VEL}[\delta ]\}}^{-1/3},\end{array}$$(7)

where r is the radius of a sphere in real space and r^* is the radius of that sphere in redshift space. Additionally, it has been found comparing the average outflow around the relevant voids with that given by the spherical expansion model that the value of g is around 0.85. This is also the value that provides the best agreement with simulations.

VEL(δ) is defined so that the peculiar velocity, V, of mass element at distance r from the centre of a spherical mass concentration (or defect) enclosing actual density contrast δ is given by

$$\begin{array}{c}\hfill V=Hr\mathrm{VEL}(\delta ),\end{array}$$(8)

where H is the Hubble constant at the time being considered. In Betancort-Rijo et al. (2006), it was shown that

$$\begin{array}{c}\hfill \mathrm{VEL}(\delta )=-\frac{1}{3}\frac{dlnD(a)}{\mathit{dlna}}\frac{\mathrm{DELK}(\delta )}{1+\delta }{\left(\frac{d}{d\delta }\mathrm{DELK}(\delta )\right)}^{-1},\end{array}$$(9)

where D(a) is the growth factor as a function of the expansion factor, a, and DELK(δ) is the inverse function of DELT(δ_l) and next (see Mo & White 1996; Sheth & Tormen 2002):

$$\begin{array}{cc}\hfill \mathrm{DELK}(\delta )& =\frac{{\delta }_c}{1.68647}(1.68647-\frac{1.35}{{(1+\delta )}^{2/3}}-\frac{1.12431}{{(1+\delta )}^{1/2}}\hfill \\ \hfill & +\frac{0.78785}{{(1+\delta )}^{0.58661}}),\hfill \end{array}$$(10)

where δ_c is the linear density contrast for spherical collapse model; for the concordance cosmology, the value at present is 1.676.

Finally, the expression of ${\displaystyle \frac{dlnD(a)}{\mathit{dlna}}}$ is:

$$\begin{array}{c}\hfill \frac{dlnD(a)}{\mathit{dlna}}\approx 1.06{\left(\frac{{(1+z)}^3}{[{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}/{\mathrm{\Omega }}_{\mathrm{m}}]}\right)}^{0.6}.\end{array}$$(11)

It is important to note that in real space, Ω_m appears in the equations only through D(a), as per Eq. (B.16)). However, we can see that in redshift space Ω_m enters the theoretical framework also through VEL(δ) (see Eq. (9)). Therefore, we would expect that the dependency of the ${\overline{n}}_v(r)$ on Ω_m is stronger in redshift than in real space.

In Figure 5, the number density of voids larger than r in redshift space is shown for the four halo simulation boxes. The agreement is still good, so we can conclude that the theoretical framework works in redshift space as well. It is important to note that the values of ${\overline{n}}_v(r)$ for the three small box catalogues (each with a different value of σ₈) are very similar for small values of r (r < 13 h⁻¹ Mpc). The differences between the three simulations begin to become important for r > 16 h⁻¹ Mpc approximately. We thus considered only these voids to constrain the cosmological parameters in the next section.

Fig. 5. Number density of voids larger than r in redshift space is shown for the theoretical framework (lines) and simulations (dots) in the left panel, while the ratio between simulations and theoretical framework is shown in the bottom panels for the Uchuu, P18, Low, and VeryLow box catalogues with a number density of $\overline{n}=3\times {10}^{-3}{h}^3$ Mpc−3. The shaded region in the bottom panel indicates delimits the region between 0.9 and 1.10 for the ratio.

In Figure 3 we have already shown the abundance of voids larger than r for SDSS and Uchuu-SDSS. The values given by the theoretical framework are shown in the same figure with a continuous black line. We can see that the agreement between the theoretical framework and Uchuu-SDSS is good for r > 14 h⁻¹ Mpc, namely, they are compatible within 1σ.

To sum up these findings, in this section we have checked that the theoretical framework successfully predicts the number density of voids larger than r for the four halo simulation boxes in real and redshift space. We have also shown that the prediction from Uchuu-SDSS voids is compatible with SDSS voids statistics within 1σ and with the prediction of the theoretical framework. Therefore, we are ready to constrain σ₈ and Ω_m and H₀ from the SDSS galaxy sample we have chosen using the theoretical framework in redshift space.

5. Bayesian analysis for cosmological parameters inference

To infer σ₈, Ω_m and H₀ from the number density of voids larger than r, we used a Bayesian analysis and Markov chain Monte Carlo (MCMC) technique to sample the posterior distribution of the considered parameter sets, Θ:

$$\begin{array}{c}\hfill \mathcal{P}(\mathrm{\Theta }|\mathcal{D})\propto \mathcal{L}(\mathcal{D}|\mathrm{\Theta })p(\mathrm{\Theta }),\end{array}$$(12)

where p(Θ) is the prior distribution, and ℒ(𝒟 ∣ Θ) is the likelihood, calculated as

$$\begin{array}{c}\hfill \mathcal{L}(\mathcal{D}|\mathrm{\Theta })={\displaystyle \prod _{i=1}^6}\frac{1}{{\sigma }_i}\exp (-\frac{{(N_{v,i}(\mathcal{D})-{\overline{N}}_{v,i}(\mathrm{\Theta }))}^2}{2{\sigma }_i^2}).\end{array}$$(13)

N_v, i(𝒟) is the number of voids within ith bin from simulations or data, ${\overline{N}}_i(\mathrm{\Theta })$ is calculated from Eq. (3) subtracting the value in the (i + 1)th bin (or r + Δr), and σ_i is the rms of N(r_i):

$$\begin{array}{c}\hfill {\sigma }_i=\mathrm{rms}({\overline{N}}_{v,i}(\mathrm{\Theta }))=\sqrt{{\overline{N}}_i(\mathrm{\Theta })/N_r},\end{array}$$(14)

where N_r represents the number of realisations for each simulation. For the simulation boxes and SDSS, N_r is equal to 1, indicating a single realisation for each set of cosmological parameters and one SDSS survey. However, in the case of Uchuu-SDSS, N_r is 32, as we have 32 lightcones of this particular type.

In this work, we use voids larger than r = 16 h⁻¹ Mpc to constrain the cosmological parameters of interest, so in Eq. (13) we consider the radius bins starting from r = 16 h⁻¹ Mpc to r = 21 h⁻¹ Mpc (i.e. 6 bins with width Δr = 1 h⁻¹ Mpc).

With Eq. (12) the best estimate of the parameters Θ may be obtained by maximising 𝒫(Θ|𝒟) with respect to the parameters Θ. To do so, we assigned to all the parameters (σ₈, Ω_m, and H₀) wide enough priors (see first column of Table 3).

To infer σ₈, Ω_m, and H₀ from SDSS and Uchuu-SDSS voids, we use MCMC sampler from Cobaya (Torrado & Lewis 2019, 2021). We assess chain convergence using the Gelman-Rubin test (Gelman & Rubin 1992). Once a tolerance of 0.01 is achieved, we consider the chains to have converged. Moreover, we discarded the first 30% of each chain as burn-in and used the mean along with the 68% two-tail equal-area confidence limit to represent the best-fit value and its uncertainty.

Additionally, we may want to combine our results with other works. In this work, we consider Planck 2018 (Planck Collaboration VI 2020) and two weak lensing results: KiDS-1000 and DESY3 (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023). To do so, we used a brand new code called CombineHarvesterFlow⁹, which allowed us to efficiently sample the joint posterior of two non-covariant experiments with a large set of nuisance parameters (Taylor et al. 2024). Specifically, CombineHarvesterFlow trains normalising flows on posterior samples to learn the marginal density of the shared parameters. Then, by weighting one chain by the density of the second flow, we can find the joint constraints.

The marginalised posteriors in the plane σ₈-H₀ (fixing Ω_m in each case to its real value; see Table 1) obtained for the four halo simulation boxes can be seen in the Appendix F. There, we show that we recover the values of σ₈ and H₀ of the simulations within 1σ (2σ) in real (redshift) space. Therefore, we can infer these cosmological parameters in Uchuu-SDSS lightcones.

6. Cosmological constraints from SDSS survey

In this section, we show the constraints obtained from SDSS redshift survey. Firstly, we show the constraints we obtain directly from the theoretical framework developed in this work and, finally, we combine our likelihood with the KIDS-1000 and DESY3 weak lensing results.

6.1. SDSS voids-only constraints

The initial contour we examine lies within the plane showcasing the core parameters of our theoretical framework: the σ₈-Γ plane, as depicted in Figure 6. In this figure, we have also included the marginalised posteriors for Uchuu-SDSS and Planck 2018.

Fig. 6. Marginalised posteriors in σ8 − Γ plane obtained for Uchuu-SDSS voids and SDSS voids voids using the maximum likelihood test with Bayesian approach with the theoretical framework developed in this work. We also show the marginalised posteriors for Planck 2018. The contours indicate the 68% (1σ) and 95% (2σ) credible intervals.

From Figure 6, it can be seen that the marginalised posterior region for Planck 2018 is entirely encapsulated within the SDSS contour at the region of 2σ, indicating that the constraints from both samples are statistically compatible in this limit. Therefore, we won’t combine SDSS with Planck 2018: the combined marginalised posteriors will be similar to the marginalised posteriors from Planck 2018.

In the left part of Figure 7 we can observe the marginalised posteriors obtained in the rest of the planes: Ω_m-H₀, σ₈-H₀, and σ₈ − Ω_m. In the right part of the same figure, we can observe the marginalised posteriors obtained from Uchuu-SDSS.

Fig. 7. Marginalised posteriors in Ωm-H0, σ8-H0, and σ8 − Ωm planes obtained for SDSS redshift survey (left) and Uchuu-SDSS lightcones (right) using the maximum likelihood test with Bayesian approach with the theoretical framework developed in this work. The contours indicate the 68% (1σ) and 95% (2σ) credible intervals.

In the first and second columns of Table 4, we can observe the constrained values directly obtained from our theoretical framework of σ₈, Ω_m, H₀, Γ, and S₈ from Uchuu-SDSS lightcones and SDSS survey, respectively. From this table, it can be seen that the uncertainties of σ₈ and Γ are much larger for SDSS than for Uchuu-SDSS because of the huge difference in the volumes between them. However, this is not the case for Ω_m and H₀. This is because σ₈ and Γ are the fundamental parameters of the theoretical framework (and S₈ depends strongly on σ₈). With these large errors, our constraints from SDSS are compatible with Planck’s within 1σ.

6.2. SDSS voids + Weak lensing

Finally, we can combine our results with KiDS-1000 and DESY3. These results are presented in Dark Energy Survey and Kilo-Degree Survey Collaboration (2023). When combining our results with other studies, such as those mentioned above, caution must be taken, as the same range for all the parameters must be considered if we want to visually compare the marginalised posteriors obtained in each case and combine the chains using CombineHarvesterFlow. This allow us to make a fair comparison of our marginalised posteriors with those of weak lensing. However, as we have seen previously, the theoretical framework used in this work depends on the parameters σ₈, Ω_m, and H₀. In weak lensing studies, A_s (DESY3) or S₈ (KiDS and KiDS-1000+DESY3) are sampled instead of σ₈, and ω_c and ω_b instead of Ω_m (KiDS and KiDS-1000+DESY3). Therefore, we have to determine which priors on σ₈ and Ω_m correspond to the aforementioned priors, ensuring compatibility across these parameters. Nevertheless, it is important to consider a significant limitation: some approximations from our theoretical framework are not valid for σ₈ values lower than 0.5, as mentioned earlier, nor for Ω_m values lower than 0.15. Thus, we have added this additional condition in the priors. This may slightly affect the process of combining chains using CombineHarvesterFlow, but, as we will see later, the effect will be very small because of the relative orientation and size of SDSS voids and weak lensing marginalised posteriors. The ranges we have used for our chains when combining with DESY3 and KiDS-1000 can be seen in the first and second columns, respectively, of Table 3. We can see that the range of H₀ used in these works is very narrow.

The marginalised posteriors in the plane σ₈ − Ω_m using these ranges can be seen in Figure 8. In this figure, the marginalised posteriors from SDSS voids with the same ranges as the different weak lensing works considered are shown, as well as each weak lensing work marginalised posteriors and the combination of SDSS voids with these three contours. We can see that SDSS voids contour is almost orthogonal to the three weak lensing contours (as expected, see Contarini et al. 2023 for more details); thus, we can anticipate that our constrained values of σ₈ and Ω_m will have smaller uncertainties than the values of the original work.

Fig. 8. KiDS-1000, DESY3, and KiDS-1000+DESY3 contours, shown in the left, middle, and right panels, respectively. The SDSS voids contours are given in the σ8 − Ωm plane, using the same parameter ranges to make a fair comparison, along with the combination of SDSS voids with each weak lensing work. In Table 3, the parameter ranges of each work are indicated.

The best-fit values of σ₈, Ω_m, H₀, and S₈ obtained from SDSS voids, KiDS-1000 and DESY3 can be seen in the first three columns of the Table 5. The combination of SDSS voids with KiDS-1000, DESY3, and KiDS-1000+DESY3 can be seen in the fourth, fifth, and sixth columns, respectively, of the same table¹⁰. We can see that the effect of combining these three weak lensing works with SDSS voids is increasing the value of σ₈, as the best-fit value of SDSS voids is very high, which implies a decrease in the best-fit value of Ω_m (from Figure 8 we can see that the combination of the three weak lensing works predicts a strong correlation between σ₈ and Ω_m, and that if one increases, the other decreases. This correlation is kept also when combining weak lensing works with SDSS voids).

We can also see from Table 5 that there is an increase in the uncertainties of H₀ when combining weak lensing with SDSS voids, which is caused by our huge uncertainties in this parameter. Specifically, the uncertainties when combining with SDSS void statistics with KiDS increase by a factor of 2–3 and a factor of 1.1 when combined with DESY3. However, the uncertainties of the rest of the parameters are decreased by a factor 2–3, approximately, with respect to the original errors of each weak lensing work.

As a final remark, we can see from Table 5 that any value of H₀ or Ω_m obtained when combining SDSS voids with the three weak lensing works are compatible with Planck 2018 within 1σ (for Ω_m, they are within 2.3σ, approximately). We can also see that S₈ value from SDSS voids + DESY3 is compatible with Planck 2018 within 1σ, but SDSS voids + KiDS-1000 and SDSS voids + KiDS-1000+DESY3 are not (noting they are compatible with Planck 2018, considering the uncertainty as 1.8σ and 1.6σ, respectively).

7. Comparison with other works about voids

In this section, we compare our constraints on σ₈, H₀, and S${}_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ with those obtained in other studies that also utilise void statistics to constrain these parameters.

The first study to compare our results with is Contarini et al. (2024), where the values of S₈ and H₀ are constrained using the void size function (the abundance of voids with radius r) predicted by the excursion set theoretical framework (Press & Schechter 1974; Sheth & Van De Weygaert 2004), employing an extension of the widely used volume-conserving model (Vdn model, Jennings et al. 2013). The constraints they derived by combining void counts with void shapes (Hamaus et al. 2020) are σ₈ = 0.809^+0.072−0.068, Ω_m = 0.308^+0.021−0.018, H₀ = 67.3^+10.0−9.1, and S${}_8=0.{813}_{-0.068}^{+0.093}$. To constrain these cosmological parameters, they used data from the BOSS DR12 redshift survey (Dawson et al. 2013); specifically, the LOWZ and CMASS target selections. The catalogues are divided into two redshift bins: 0.2 < z ≤ 0.45 and 0.45 < z < 0.65. This sample encompasses a physical volume approximately 60 times larger than the one used in this work. According to Reid et al. (2015), the CMASS volume is 5.1 Gpc³, and the LOWZ volume is 2.3 Gpc³, resulting in a total volume of 2.3 × 10⁹ h⁻³ Mpc³. If we consider the scaling of errors with survey volume (as demonstrated in Appendix G) and apply this scaling to a redshift survey like BOSS DR12 using our theoretical framework, we can estimate how much the uncertainties would decrease. For Γ (a parameter that was not constrained in Contarini et al. 2024), our errors would scale by a factor of 60^−1/2.5 = 0.194, making them five times smaller than those obtained in this work with the SDSS sample volume. For σ₈, the uncertainties would scale by 60^−1/7 = 0.557, reducing them by a factor of 1.79. This would result in uncertainties approximately 30% larger than those reported in Contarini et al. (2024). For S₈, the uncertainties would scale by 60^−1/6 = 0.505, reducing them by a factor of 1.98. This again implies uncertainties roughly 30% larger than those in Contarini et al. (2024). However, for H₀ and especially Ω_m, there would be no significant improvement due to the limited scaling of these parameters with survey volume. As a result, the uncertainties for H₀ and Ω_m would remain significantly larger than those obtained in Contarini et al. (2024). It is important to note that the constraints provided in the aforementioned study are derived from the combination of the Vdn model and void shapes, whereas the constraints presented in this work are based solely on void statistics.

Another study to compare our results with is Sahlén et al. (2016), where galaxy cluster and void abundances are combined using extreme-value statistics. This approach focuses on a large galaxy cluster and a void aligned with the Cold Spot (CS) in the CMB (referred to as the CS void) (Finelli et al. 2016). Using this method, they obtained a constraint of σ₈ = 0.95 ± 0.21 for a flat ΛCDM universe. This result is consistent with values derived from the CMB as well as weak lensing studies.

We now compare our results with those from CMB, weak lensing, Cepheid-based measurements, and tip of the red giant branch experiments. Figures 9 and 10 display the constraints on σ₈, H₀, and S₈ obtained in this work alongside those from the aforementioned studies. As discussed in the introduction, there is a significant discrepancy between the constraints on S₈ (and σ₈) and H₀ derived from high-redshift and low-redshift cosmological probes. This discrepancy is evident in Figures 9 and 10, where it is clear that the constraints from CMB and weak lensing are not consistent with those from Cepheids and SNIa measurements.

Fig. 9. Comparison between recent constraints on the parameters σ8 (left panel) and H0 (right panel) from different cosmological probes. The error bars represent 68% confidence intervals. The black error bars represent the values constrained in this work. The references of the rest of the works, from top to bottom in left panel are: Dark Energy Survey and Kilo-Degree Survey Collaboration (2023), Heymans et al. (2021), and Planck Collaboration VI (2020), and for the right panel: Freedman et al. (2020), Riess et al. (2022), Brout et al. (2022a,b), and Planck Collaboration VI (2020).

Fig. 10. Comparison between recent constraints on S8 from different cosmological probes. The error bars represent 68% confidence intervals. The black error bars represent the values constrained in this work. The references of the rest of the works, from top to bottom are: Dark Energy Survey and Kilo-Degree Survey Collaboration (2023), Li et al. (2023), and Planck Collaboration VI (2020).

In the σ₈ panel of Figure 9, we observe that combining our results with KiDS-1000 or KiDS-1000+DESY3 yields σ₈ values that are compatible with Planck 2018 within 1σ. However, combining our results with DESY3 alone produces a combined σ₈ value that is incompatible with both Planck 2018 and DESY3 within 1σ. Additionally, the combined H₀ values from all three weak lensing studies are consistent within 1σ with the H₀ measurements obtained from Cepheids+SNIa.

In Figure 10, we can see a comparison of the constrained S₈ values obtained by combining SDSS voids with KiDS-1000 or KiDS-1000+DESY3 with CMB and weak lensing constraints. We can see that the constraints from SDSS combined with the three weak lensing studies considered in this work are consistent within 1σ. When we combine SDSS voids with DESY3, the resulting S₈ value is also compatible within 1σ with Planck 2018, but the rest of the combinations are not.

To conclude this section, we estimate the potential reduction in errors, considering only void statistics, with future redshift surveys such as Euclid (Euclid Collaboration 2020) and the Dark Energy Spectroscopic Instrument (DESI) experiment (DESI Collaboration 2016a,b). Euclid aims to investigate the Universe’s expansion history and the evolution of large-scale structures by measuring the shapes and redshifts of galaxies, covering 15 000 deg² of the sky and extending to redshifts of approximately z = 2.

Similarly, DESI will perform the Bright Galaxy Survey (BGS), targeting over 10 million galaxies in the redshift range 0 < z < 0.6, as well as a dark-time survey of 20 million luminous red galaxies (LRGs), emission-line galaxies (ELGs), and quasars (Hahn et al. 2023). With a projected footprint of 14 000 deg² and a broader redshift range, DESI is expected to achieve a precision 1–2 orders of magnitude higher than that of existing surveys, such as SDSS and BOSS.

We can estimate the reduction in uncertainties for the main parameters of the theoretical framework (σ₈ and Γ) when using DESI Y1 by calculating the ratio between the volumes of the SDSS and DESI surveys to the power of a constant value which has been calculated for each parameter using the constraints from a different number of lightcones (see Appendix G for more details). For σ₈, the uncertainties would decrease as

$$\begin{array}{c}\hfill \mathrm{\Delta }\approx {\left(\frac{V_{\mathrm{BGS},\mathrm{Y}1}}{V_{\mathrm{SDSS}}}\right)}^{-1/7}\sim {\left(\frac{1.2\times {10}^9{h}^{-3}{\mathrm{Mpc}}^3}{40\times {10}^6{h}^{-3}{\mathrm{Mpc}}^3}\right)}^{-1/7}\sim 0.62;\end{array}$$(15)

in other words, the uncertainties for these two parameters would be 1.63 times smaller than the uncertainties obtained with the sample of SDSS used in this work.

For Γ the decrease would be

$$\begin{array}{c}\hfill \mathrm{\Delta }\approx {\left(\frac{V_{\mathrm{BGS},\mathrm{Y}1}}{V_{\mathrm{SDSS}}}\right)}^{-1/2.5}\sim {\left(\frac{1.2\times {10}^9{h}^{-3}{\mathrm{Mpc}}^3}{40\times {10}^6{h}^{-3}{\mathrm{Mpc}}^3}\right)}^{-1/2.5}\sim 0.26,\end{array}$$(16)

which implies that the uncertainties for this parameter would be 3.9 times smaller than the uncertainties obtained with the sample of SDSS used in this work.

Therefore, the uncertainties for σ₈ (and S₈) obtained with the Bright Galaxy Survey from DESI Y1 will be approximately a factor 1.7 lower than those obtained in this work with the SDSS sample used in this work, and a factor 3.9 for Γ. These uncertainties will be even lower when the full DESI survey is complete.

8. Summary and future developments

In this work, we have made use of the theoretical framework developed in Betancort-Rijo et al. (2009) and re-calibrated the expression for the number density of voids greater than r (see Eq. (3)). We used the Uchuu halo simulation box with a number density of halos equal to 3 × 10⁻³ h³ Mpc⁻³, which has a greater volume than the simulation used in the aforementioned work (V = 2000³ h⁻³ Mpc³). The most important results obtained in this work are as follows:

• We have proven a compatibility within 1σ among the number density of voids greater than r of SDSS galaxies with M_r < −20.5 (where M_r is the absolute magnitude in r-band), z_min = 0.02, and z_max = 0.132, as well as Uchuu-SDSS galaxies with the same characteristics.

• We have demonstrated that our theoretical framework can successfully predict the abundance of large voids for the four halo simulation boxes with different values of σ₈ (σ₈ = {0.8159, 0.8102, 0.75, 0.65}) and Uchuu-SDSS lightcones used in this work,

• We used a Bayesian analysis for all halo simulation boxes and calculated the contours for 68% (1σ) and 95% (2σ) credible intervals in σ₈-H₀ plane (fixing Ω_m in each case to its real value; see Table 1). We proved that we can recover the real values of σ₈ and H₀ of each halo simulation box within 1σ in real space and 2σ in redshift space,

• We checked that using MCMC sampler from Cobaya, we successfully recover the values of these parameters for Uchuu-SDSS lightcones. The recovered values are ${\sigma }_8=1.{296}_{-0.592}^{+0.621}$, ${\mathrm{\Omega }}_{\mathrm{m}}=0.{419}_{-0.152}^{+0.172}$, H${}_0=72.{35}_{-11.85}^{+12.24}$, $\mathrm{\Gamma }=0.{250}_{-0.088}^{+0.096}$, and S${}_8=1.{61}_{-0.93}^{+1.05}$.

• We calculated the contours in σ₈ − Ω_m, σ₈-H₀, and H₀ − Ω_m planes for the same ranges as mentioned above for SDSS void statistics. The contours obtained for SDSS are much broader than the ones obtained for halo simulation boxes (and Uchuu-SDSS lightcones) because of the huge difference of volume (the SDSS volume is approximately 25 times smaller than the three small simulation boxes and 200 times smaller than Uchuu box volume), which means that there are far fewer voids in SDSS than in the boxes. The constrained values obtained from SDSS are ${\sigma }_8=1.{520}_{-0.441}^{+0.416}$, ${\mathrm{\Omega }}_{\mathrm{m}}=0.{459}_{-0.184}^{+0.184}$, H${}_0=71.63{\pm }_{-11.77}^{+12.60}$, $\mathrm{\Gamma }=0.{270}_{-0.100}^{+0.094}$, and S₈ = 1.87${}_{-0.76}^{+0.59}$. It is important to remark that these constrained values have been obtained supposing that the ratio of Ω_b/Ω_m is constant and given by Planck 2018.

• Next, we combined our SDSS voids constraints with KiDS-1000, DESY3, and KiDS-1000+DESY3 (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023). The results obtained when combining SDSS voids with KiDS-1000+DESY3 are: ${\sigma }_8=0.{858}_{-0.040}^{+0.040}$, ${\mathrm{\Omega }}_{\mathrm{m}}=0.257{\pm }_{-0.020}^{+0.023}$, H${}_0=74.{17}_{-4.66}^{+4.66}$, and S₈ = 0.794${}_{-0.016}^{+0.016}$. We see that when we combine our results with weak lensing works, we significantly improve the uncertainties of all the parameters (with the exception of H₀) with respect to the uncertainties of only KiDs-1000+DESY3. This is because the marginalised posteriors in the σ₈ − Ω_m plane obtained in this work and those from KiDS-1000 are almost orthogonal.

• We have seen that any value of H₀ or Ω_m obtained when combining SDSS voids with the three weak lensing works are compatible with Planck 2018 within 1σ (for Ω_m, they are within approximately 2σ). We have also seen that S₈ value from SDSS voids + DESY3 is compatible with Planck within 1σ, but SDSS voids + KiDS-1000 and SDSS voids + KiDS-1000 + DESY3 are not (the former is compatible within 2.7σ and the latter within 2σ).

• Finally, we compared our results with those obtained in other works where voids have also been used to constrain the same cosmological parameters constrained in the current work, along with parameters constrained via CMB, galaxy clustering, weak lensing, and Type Ia supernova measures. We have checked that when our sample is combined with those of weak lensing studies, we obtained slightly smaller uncertainties than in other works focussed solely on voids. However, when we do not combine our results with any work, our uncertainties are very large because the volume of the samples of Uchuu-SDSS and SDSS used in this work is relatively small, in comparison to the volumes of redshift surveys used in other works of voids. For example, in Contarini et al. (2024), the BOSS redshift survey was used. If we used this redshift survey, we would have predicted that our uncertainties (without combining with any work) for σ₈ and S₈ would be approximately a 30% larger than the ones given in Contarini et al. (2024) (where the authors combined void statistics with void shapes).

• Therefore, as we have demonstrated in Section 7, the most important limitation we face when constraining the cosmological parameters considered in this study is the small volume of the redshift survey sample we use, which is approximately ∼ 40 × 10⁶ h⁻³ Mpc³. When we use a survey such as DESI Y1, we would have estimated uncertainties for Γ (σ₈) that are approximately 4 (1.5) times smaller than the uncertainties obtained with the sample of SDSS used in this work. If we combined our study with other, additional works, then the decrease in the errors would be greater, as a result of increasing the volume of the sample of galaxies used.

Data availability

The Uchuu halo and galaxy boxes and the 4 box catalogs at redshift z = 0.092, as well as the 32 Uchuu-SDSS galaxy lightcones, the SDSS catalogue, and the void catalogues from all the previous galaxy and halo catalogues used in this work are available at: http://www.skiesanduniverses.org/Simulations/Uchuu/. This link includes information on how to read the data and column description. For a list and brief description of the available halo, Uchuu-SDSS and SDSS void catalogues columns, see Appendix I.

Acknowledgments

E. Fernández-García acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033. EFG, FP and AK thanks support from the Spanish MICINN PID2021-126086NB-I00 funding grant. T.I has been supported by IAAR Research Support Program in Chiba University Japan, MEXT/JSPS KAKENHI (Grant Number JP19KK0344, JP21H01122, and JP23H04002), MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (JPMXP1020200109 and JPMXP1020230406), and JICFuS. We thank Instituto de Astrofisica de Andalucia (IAA-CSIC), Centro de Supercomputacion de Galicia (CESGA) and the Spanish academic and research network (RedIRIS) in Spain for hosting Uchuu DR1, DR2 and DR3 in the Skies & Universes site for cosmological simulations. The Uchuu simulations were carried out on Aterui II supercomputer at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan, and the K computer at the RIKEN Advanced Institute for Computational Science. The Uchuu Data Releases efforts have made use of the skun@IAA_RedIRIS and skun6@IAA computer facilities managed by the IAA-CSIC in Spain (MICINN EU-Feder grant EQC2018-004366-P). The other cosmological simulations were carried out on the supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ID: hp220173, hp230173, and hp230204). We also thank Joe Zuntz and Anna Porredon for their feedback on our cosmological results. Finally, we are thankful to Peter Taylor for instructing us on the the use of CombineHarvesterFlow, for verifying our results and for giving us feedback on the paper.

References

Appendix A: Void finder code description

In this appendix, we outline the custom C++ code developed to identify maximal non-overlapping spheres (voids) from the candidates provided by DIVE. The code utilizes OpenMP for parallelization. A summary of the implementation and its workflow is illustrated in Figure A.1.

Fig. A.1. Void identification process.

The code begins by including necessary headers and defining a structure for spheres with center coordinates (X, Y, Z) and radius R. It also includes a utility function to check overlaps: two spheres, s₁ and s₂, overlap if the distance between their centers satisfies:

$$\begin{array}{c}\hfill {(X_1-X_2)}^2+{(Y_1-Y_2)}^2+{(Z_1-Z_2)}^2<{(R_1+R_2)}^2\end{array}$$(A.1)

The main function performs the following steps:

1. Reads the input file containing DIVE candidates.

2. Sorts candidates by descending radius, R.

3. Initializes the void list with the largest candidate.

4. terates through the remaining candidates, checking overlaps using Eq. (A.1):

   • If no overlap is detected, the candidate is retained as a void.

   • If overlap is detected, the candidate is discarded.

5. Outputs the identified voids to a file.

The resulting output provides the centers and radii of all voids found in the halo or galaxy sample.

Appendix B: Void statistics theoretical framework in detail

In this appendix we explicitly write all necessary quantities in order to calculate P₀(r) (VPF) and ${\overline{n}}_v(r)$. This theoretical framework can also be consulted in Betancort-Rijo et al. (2009). The first step to calculate all the terms that appear in Eq. (3) is to develop Eq. (6). In Patiri et al. (2006a), it is shown that for dark matter halos, u can be written as

$$\begin{array}{c}\hfill u=[\overline{n}V(1+\delta )][1+{\delta }_{\mathit{ns}}],\end{array}$$(B.1)

where $\overline{n}$ denotes the mean number density of those halos in the sample (usually halos larger than some given mass), V is the volume of the sphere and δ the actual enclosed density contrast within the sphere. The first factor on the right-hand side of the equation is the integral of the probability density within the sphere for halos tracing the mass (i.e. no bias, which is true in the very low mass limit). In general, halos are biased tracers of the underlying mass distribution, due to the initial clustering of the protohalos before they move along with mass (i.e. statistical clustering). The second term of the equation accounts for this biasing. In Patiri et al. (2006a) an approximation for this bias is obtained as a function of the linear enclosed density contrast within the sphere (δ_l):

$$\begin{array}{c}\hfill 1+{\delta }_{\mathit{ns}}({\delta }_l)=A(m)e^{-b(m){\delta }_l^2}\forall {\delta }_l\le -1,\end{array}$$(B.2)

where A(m), b(m) are coefficients mainly depending on the halo mass. In Patiri et al. (2006a) is demonstrated that, using Zeldovich approximation, P_n(r) can be rewritten as

$$\begin{array}{c}\hfill P_n(r)={\displaystyle {\int }_{-\infty }^{1.6}}P({\delta }_l,r)\frac{{[u({\delta }_l)]}^n}{n!}{e}^{[-u({\delta }_l)]}d{\delta }_l,\end{array}$$(B.3)

where u(δ_l) is now a function of δ_l through the dependence of the actual density contrast, δ, on its linear counterpart, δ_l:

$$\begin{array}{c}\hfill u({\delta }_l)=(\overline{n}V[1+\delta ({\delta }_l,r)])[1+{\delta }_{\mathit{ns}}({\delta }_l)],\end{array}$$(B.4)

where δ(δ_l, r) is basically the relationship between the actual and the linear density contrast within a sphere as given by the standard spherical collapse model, except for a small correcting term depending on r.

Doing some manipulation, P_n(r) can be rewritten as

$$\begin{array}{c}\hfill P_n(r)=\frac{1}{n!}{\displaystyle {\int }_{-7}^0}P({\delta }_l,r){[u({\delta }_l)]}^n{e}^{[-u({\delta }_l)]}d{\delta }_l,\end{array}$$(B.5)

where u(δ_l) is now

$$\begin{array}{c}\hfill u({\delta }_l)=\overline{n}V[1+DELF({\delta }_l,r)]A{e}^{-b{\delta }_l^2}.\end{array}$$(B.6)

DELF is a function of δ_l and r that gives the mean actual density contrast within a sphere with radius r with enclosed linear density contrast δ_l:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill DELF({\delta }_l,r)& =\hfill \\ \hfill \frac{1+DELT({\delta }_l)}{|1-(4/21){[1+DELT({\delta }_l)]}^{2/3}{[\sigma (r{[1+DELT({\delta }_l)]}^{1/3})]}^2|}\end{array},\end{array}$$(B.7)

where DELT(δ_l) denotes the relationship between the actual and linear enclosed density contrasts in the spherical collapse model (Patiri et al. 2006a):

$$\begin{array}{c}\hfill 1+DELT({\delta }_l)\approx {(1-0.607{\delta }_l)}^{-1.66},\end{array}$$(B.8)

where σ(Q) is the rms of the linear density contrast on a sphere with Lagrangian radius Q. In this equation, σ(Q) is evaluated at Q = r[1 + DELT(δ_l)]^1/3. A, b in Eq. (B.6) are also functions of δ_l (see Rubiño-Martín et al. 2008, for more details) given by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A& =A(m,Q=r{[1+DELT({\delta }_l)]}^{1/3})\times \hfill \\ \hfill & \times {\left[\frac{D(z){\sigma }_8}{0.9}\right]}^{0.88}{\left(\frac{\mathrm{\Gamma }}{0.21}\right)}^{0.174}\hfill \end{array},\end{array}$$(B.9)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b& =b(m,Q=r{[1+DELT({\delta }_l)]}^{1/3})\times \hfill \\ \hfill & \times {\left[\frac{D(z){\sigma }_8}{0.9}\right]}^{-2.55}{\left(\frac{\mathrm{\Gamma }}{0.21}\right)}^{-0.69}\hfill \end{array},\end{array}$$(B.10)

where A(m, Q), b(m, Q) are functions of the mass of the objects and the Lagragian radius of the regions being considered:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A(m,Q)& =[1.577-0.298\left(\frac{Q}{8}\right)]-\hfill \\ \hfill & -[0.0557+0.0447\left(\frac{Q}{8}\right)]lnm-\hfill \\ \hfill & -[0.00565+0.0018\left(\frac{Q}{8}\right)]{[lnm]}^2\hfill \end{array},\end{array}$$(B.11)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b(m,Q)& =[0.0025-0.00146\left(\frac{Q}{8}\right)]+\hfill \\ \hfill & +[0.121-0.0156\left(\frac{Q}{8}\right)]m^{0.335+0.019Q/8}\hfill \end{array},\end{array}$$(B.12)

$$\begin{array}{c}\hfill m=\frac{1.05M}{3.4866\times {10}^{11}{h}^{-1}{M}_{\odot }}\left(\frac{0.3}{{\mathrm{\Omega }}_{\mathrm{m}}}\right),\end{array}$$(B.13)

where D(z) is the linear growth factor normalised to be 1 at present and M is the minimum mass of halos when the matter density field is traced by the halo positions (or central halos). However, the matter density field of Uchuu is traced by the positions of its sub-structures (or subhalos), so the mass to be put in the theoretical framework is not m, but m_g (see Appendix C.3 for a detailed explanation of how one can calculate m_g).

We already have all terms in order to calculate u(δ_l) from Eq. (B.5), so all we need now is to calculate P(δ_l, r), which is the probability function of the linear density contrast within an Eulerian space, and can be calculated as (Betancort-Rijo & López-Corredoira 2002):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P({\delta }_l,r)& =\frac{exp[(-1/2){\delta }_l^2/{(\sigma (r{[1+DELF({\delta }_l,r)]}^{1/3}))}^2]}{\sqrt{2\pi }}\times \hfill \\ \hfill & \times {[1+DELF({\delta }_l,r)]}^{-[1-(\alpha /2)]}\times \hfill \\ \hfill & \times \frac{d}{d{\delta }_l}\left(\frac{{\delta }_l}{\sigma (r{[1+DELF({\delta }_l,r)]}^{1/3})}\right)\hfill \end{array},\end{array}$$(B.14)

where

$$\begin{array}{c}\hfill \alpha ({\delta }_l,r)=0.54+0.173ln\left(\frac{r{[1+DELF({\delta }_l,r)]}^{1/3})}{10}\right),\end{array}$$(B.15)

and

$$\begin{array}{cc}\hfill \sigma (Q)& =\sigma (Q,\mathrm{\Gamma })\approx {\sigma }_{8}A(\mathrm{\Gamma })Q^{-B(\mathrm{\Gamma })-C(\mathrm{\Gamma })Q},\hfill \end{array}$$(B.16)

$$\begin{array}{cc}\hfill A(\mathrm{\Gamma })& =2.01+3.9\mathrm{\Gamma },\hfill \end{array}$$(B.17)

$$\begin{array}{cc}\hfill B(\mathrm{\Gamma })& =0.2206+0.361{\mathrm{\Gamma }}^{1.5},\hfill \end{array}$$(B.18)

$$\begin{array}{cc}\hfill C(\mathrm{\Gamma })& =0.182+0.0411ln\mathrm{\Gamma }.\hfill \end{array}$$(B.19)

This fit is valid for Q ≥ 3h⁻¹Mpc and 0.1 ≥ Γ ≥ 0.5.

It is important to remark that Eq. (B.16) is only valid for z = 0. If we want to calculate this function in a different redshfit, then we have to make the following change:

$$\begin{array}{c}\hfill {\sigma }_8\to {\sigma }_8\frac{D(z)}{D(z=0)},\end{array}$$(B.20)

where D(z) is the linear growth factor of density fluctuations in the model under consideration.

Finally, it’s important to write the equation for the root mean square (rms) of the error of the estimation of the VPF. The expression we use is the following:

$$\begin{array}{c}\hfill rms{(P_0(r))}^2=\frac{9.2[1-\omega {\overline{n}}_v(r)]P_0{(r)}^2}{N(r)}+\frac{P_0(r)}{N_{\mathit{spheres}}},\end{array}$$(B.21)

where N_spheres is the number of random spheres used in order to estimate P₀(r) for simulations, N(r) is the number of voids larger than r. This equation is a modification of the expression given in Betancort-Rijo (1992) and Patiri et al. (2006b). We have added a second term, P₀(r)/N_spheres to take into account the error due to the finite number of trial random spheres used to calculate the VPF. The first term of Eq. (B.21) takes into account the finite volume of the sample.

Then, ω in Eq. (B.21) is given by

$$\begin{array}{cc}\hfill \omega =\frac{32\pi }{3}{\overline{R}}^3& [1+2.73{\overline{n}}_v(r)\frac{32\pi }{3}{\overline{R}}^3\times \hfill \\ \hfill & \times {(1-\frac{3}{4}\frac{{\overline{n}}_v{(r)}^{-1/3}}{2\overline{R}}+\frac{1}{8}{\left(\frac{{\overline{n}}_v{(r)}^{1/3}}{2\overline{R}}\right)}^2)]}^{-1},\hfill \end{array}$$(B.22)

where $\overline{R}$ is the mean radius of all voids larger than the minimum voids considered.

Appendix C: Calculation of m_g

As shown in Appendix B, P_n(r) (in concrete, u), depends on m (M normalised according to Eq. (B.13)), which is the minimum mass of halos when the matter density field is traced by the halo positions. However, we want to use this theoretical framework to infer cosmological parameters from Uchuu, which matter density field is traced by its sub-structures (subhalos), and from SDSS galaxies, which are also substructures of halos. In this case, the mass to be used in the theoretical framework is not m but m_g (i.e. the minimum mass of subhalos –including halos).

These two masses, M and M_g, can be calculated directly from the simulations used. For example, M can be calculated imposing $\overline{n}(>M)={\overline{n}}_{\mathrm{sample}}$, where ${\overline{n}}_{\mathrm{sample}}$ is the number density of the sample considered. In our case, ${\overline{n}}_{\mathrm{sample}}$ = 3×10⁻³h³Mpc⁻³. Note that in this case $\overline{n}(>M)$ is the cumulative halo mass function. M_g can be calculated the same way, but taking also into account the subhalos in $\overline{n}(>M_g)$.

Therefore, in a simulation such as Uchuu, with halos and subhalos, we can calculate both M and M_g and establish a relationship between both. This was already done in Betancort-Rijo et al. (2009) (see Eq. (16) from that article). However, in this work we have recalibrated this expression and found that m_g can be approximately related to m through the following equation:

$$\begin{array}{c}\hfill m_g=1.{058}^{\sigma (m)}m,\end{array}$$(C.1)

where σ(m) is the rms linear density fluctuation on scale m. It can be shown that, in the interval of masses and cosmological parameters considered in this work, σ(m) depends on m, Γ and σ₈ the following way:

$$\begin{array}{c}\hfill \sigma (m)=2.035{\left(\frac{m}{4.039}\right)}^{-\frac{0.475}{3}}{\left(\frac{\mathrm{\Gamma }}{0.1754}\right)}^{0.27}\left(\frac{{\sigma }_8}{0.8159}\right),\end{array}$$(C.2)

Additionally, an expression for m_g depending on the ${\overline{n}}_{\mathit{sample}}$, σ₈ and Γ can be found:

$$\begin{array}{c}\hfill m_g=4.74{\left(\frac{{\sigma }_8}{0.8159}\right)}^{0.402}{\left(\frac{\mathrm{\Gamma }}{0.1754}\right)}^{0.109}\left(\frac{F({\overline{n}}_{\mathit{sample}})}{4.047}\right),\end{array}$$(C.3)

where

$$\begin{array}{c}\hfill F(\overline{n})=-\frac{1.267\times {10}^{-8}}{{\overline{n}}^2}+\frac{1.289\times {10}^{-2}}{\overline{n}}-0.248.\end{array}$$(C.4)

This last function, $F(\overline{n})$ defines the dependence of the mass m with the number density of the sample, $\overline{n}$. It is important to remark that ${\overline{n}}_{\mathit{sample}}$ is the number density of galaxies corrected from completeness in the case we use redshift surveys or lightcones.

From Eq. (C.3), we can see that we obtain m_g = 4.74 for Uchuu, which is the mass obtained using directly the simulations (for n_g = 3 × 10⁻³h³Mpc⁻³, we take the the N = n_gV most massive halos, where V is the Uchuu volume, and fix m_g to the minimum mass of this N halos) using the normalisation given by Eq. (B.13). This equation is a good approximation to m_g as a function of cosmological parameters and ${\overline{n}}_{\mathit{sample}}$, while the abundance matching provides it exact value. However, the dependence on σ₈ and Γ is quite small and it is swamped by small error in expressions (B.9-B.12) (which have been fitted to the results of a complex procedure) that have a much stronger dependence on σ₈ and Γ. Thus, we have found that the best agreement between the theoretical framework and the simulations are obtained holding m_g fixed for given values of ${\overline{n}}_{\mathit{sample}}$.

The value of m_g for the four halo boxes has been calculated using the chi square test with the VPF defined in real space¹¹ for r > 14h⁻¹Mpc, r < 21h⁻¹Mpc and Δr = 1h⁻¹Mpc (i.e 8 bins). We have minimised the sum of ${\chi }^2={\chi }_{\mathrm{Uchuu}}^2+{\chi }_{\mathrm{P}18}^2+{\chi }_{\mathrm{P}18\mathrm{LowS}8}^2+{\chi }_{\mathrm{P}18\mathrm{VeryLowS}8}^2$, and have obtained an optimal value of m_g = 4.90 with χ² = 41.67 and ${\chi }_{\nu }^2={\chi }^2/\nu$ = 1.49 (where ν is the degree of freedom, ν = 4 × 7 in our case). Therefore, m_g = 4.90 is the value we will use in this work.

Appendix D: Validation of the theoretical framework in real and redshift space

As already mentioned, the theoretical framework predicts the void probability function (VPF) in real and redshift space from first principles and establish a relationship between the VPF and the number density of voids larger than r. This relationship between the VPF and the number density of voids larger than r does not arise from first principles but is an assumption that, although not exact, is quite accurate (see Appendix E for more details).

In Section 4 and Appendix B, we give all the expressions needed to compute the VPF and m_g, although in Appendix C we calculate an universal mass for the 4 simulation boxes used in this work as the dependence of m_g with σ₈ is small.

Therefore, we are in the position of checking if the VPF can reproduce the VPF from the 4 simulation boxes in real and redshift space. This can be seen in Figure D.1 (real space) and D.2 (redshift space). We can see that the values given by all the simulations are compatible within 1σ with the values predicted by the theoretical framework, both in real and redshift space. Therefore, we can conclude that the theoretical framework, which predicts the VPF based on first principles, is valid.

Fig. D.1. VPF in real space is shown for the theoretical framework (lines) and simulations (dots) given in the top panel, while the ratio between simulations and theoretical framework is shown in bottom panels for the Uchuu, P18, Low, and VeryLow box catalogues with number density $\overline{n}=3\times {10}^{-3}{h}^3$Mpc−3. Shaded region in bottom panel indicates delimits the region between 0.9 and 1.10 for the ratio.

Fig. D.2. VPF in redshift space is shown for the theoretical framework (lines) and simulations (dots) given in the top panel, while the ratio between simulations and theoretical framework is shown in bottom panels for the Uchuu, P18, Low and VeryLow box catalogues with number density $\overline{n}=3\times {10}^{-3}{h}^3$Mpc−3. Shaded region in bottom panel indicates delimits the region between 0.9 and 1.10 for the ratio.

Appendix E: Calibration of the relationship between the VPF and the number density of voids larger than r

The theoretical framework presented in this work is based in the model presented in Betancort-Rijo et al. (2009). This model was calibrated with Millennium Run simulations (Springel et al. 2005), which follows the evolution of 1010 dark matter particles in a periodic box of 500 h⁻¹ Mpc on a side with a mass resolution per particle of 8.6×10⁸ h⁻¹ M_⊙. This volume is extremely small in comparison with Uchuu’s volume (Uchuu’s volume is 64 times larger than Millennium Run’s), so the void statistics calculated with this box weren’t as accurate as the ones obtained with Uchuu, which translates in an imprecision in the calibration of the theoretical framework.

The parameters of the theoretical framework that were calibrated with simulations were μ, α and β, defined by

$$\begin{array}{c}\hfill {\overline{n}}_v(r)=\frac{\mu \mathcal{K}(r)}{V}{e}^{-\alpha \mathcal{K}(r)[1-\beta K(r)]},\end{array}$$(E.1)

along with the mass, m_g, to be used in all the halo and galaxy boxes, as well as lightcones and SDSS. The calculation of m_g is explained in Appendix C.3.

In order to find the values of the coefficients μ, α, and β, we have used the chi square test with the number density of voids larger than r predicted by the theoretical framework used in this work, letting these coefficients be free parameters, and the values given by the simulations. We have minimised sum of the chi square values of each halo box: ${\chi }^2={\chi }_{\mathrm{Uchuu}}^2+{\chi }_{\mathrm{P}18}^2+{\chi }_{\mathrm{P}18\mathrm{LowS}8}^2+{\chi }_{\mathrm{P}18\mathrm{VeryLowS}8}^2$ and we have obtained the optimal value of μ = 0.603, α = 1.792 and β= 0 with χ² = 26.99 and ${\chi }_{\nu }^2={\chi }^2/\nu$ = 1.12 (where ν is the degree of freedom, ν = 4 × 6 in our case). Therefore, these are the values of μ, α, and β that we use in this work.

Appendix F: Constraints of cosmological parameters using voids statistics in the distribution of halos and galaxies in mocks with theoretical framework

In this appendix, we show the marginalised posteriors obtained with the theoretical framework presented in Appendix B using the maximum likelihood test with Bayesian approach. We compare it with the results of the halo simulation boxes. We do this in real and redshift space.

F.1. Real Space

We can see the marginalised posteriors obtained in the for halo boxes in real space in the first column of Figure F.1. In this figure, we can see that all fiducial values (black points) are inside the 1σ contour.

Fig. F.1. Constraints from Nv, i identified in Uchuu (first row), P18 (second row), Low (third row) and VeryLow (row) halo simulation boxes in real (first column) and redshift (second column) spaces using Maximum Likelihood test. The contours indicate the 68% (1σ) and 95% (2σ) credible intervals. The black dots are the real values of σ8 and H0 of each simulation.

Therefore, we successfully recover the values of σ₈ and Γ of the simulations with the theoretical framework, and we can take one more step to do this study in redshift space.

F.2. Redshift space

The marginalised posteriors obtained with maximum likelihood test for the halo simulation boxes in redshift space can be seen in the second column of Figure F.1. In that figure, we can see that we obtain a small contour for Uchuu in redshift space, unlike in real space. Additionally, we recover Planck 2018 values of σ₈ and H₀ within 1σ for the P18 box catalogue, and within 2σ for Low and VeryLow. Therefore, in redshift space we have not successfully recovered the real values of the parameters in simulations within 1σ as we have in real space. This can be due to the transformation we have used in order to transform the theoretical framework from real to redshift space (see Eq. (7)). However, we can recover the real values within 2σ for the four halo simulation boxes.

Appendix G: Scaling of errors with volume

An important aspect to consider is the scaling of errors with the sample volume, specifically: how the errors decrease as the sample volume increases.

We can observe the marginalised posterior in the σ₈ − Ω_m plane in Figure G.1. From this figure, it is clear that the marginalised posteriors (and, consequently, the constrained values of each parameter) significantly decrease when a larger sample volume is considered. Additionally, we see that the marginalised posterior for a single random Uchuu-SDSS lightcone void is similar to that obtained for SDSS voids.

Fig. G.1. VPF (left panel) and number density of voids larger than r (right panel) predicted by the theoretical framework developed in this work (continuous lines) and measured in Uchuu galaxy box (points) for different galaxy number densities, ng, at redshift z = 0.092.

In Table G.1, we present the constrained values for the voids from a single Uchuu-SDSS lightcone (second column), 5 lightcones (third column), 10 lightcones (fourth column), and 32 lightcones (fifth column) for the main parameters of the theoretical framework: σ₈ and Γ. The last column shows the power that best fits the relationship y = x^1/n, where x is the volume ratio between each column (all possible combinations, V₂/V₁ with V₁ < V₂) and y is the ratio of uncertainties. We performed two different fits: one for the upper limit uncertainties and another for the lower limit uncertainties. The value of n was averaged from the two fits and then rounded to the nearest integer. However, for Γ, we obtained n = −2.5, so no rounding was applied in this case. The fits were conducted using non-linear least squares.

Therefore, we can observe that the scaling of errors with volume differs for each parameter. The largest scaling is found for Γ, where the volume scaling is V^−1/2.5. This means that if we have two samples and one has twice the volume of the other, the errors in the larger sample would be 0.758 times the errors in the smaller sample (i.e. a factor of 1.32 smaller). However, the scaling for σ₈ (and S₈) is less favourable: V^−1/7. Using the same example, if one sample has twice the volume of the other, the errors in the larger sample would be 0.906 times the errors in the smaller sample (i.e. a factor of 1.104 smaller).

One way to check the robustness of the values of n is by recalculating the inferences for a different random group of 10 lightcones and comparing the resulting n values. Instead of performing the fit with all cases from Table G.1, we consider one lightcone and one group of 10 lightcones, then another lightcone and the other group of 10 lightcones, and finally compare the n values obtained in each case. In the first two columns of Table G.2, we show the inferences obtained for each group (A and B) and the corresponding n values (n_A and n_B). The error in n, denoted as Δn, is calculated as the absolute difference between n_A and n_B.

Therefore, we can conclude that n_σ8 ∼ −7 ± 2 and n_Γ ∼ −2.5 ± 0.1, so the calculation of n, specially for Γ, is quite robust.

Appendix H: Theoretical framework with different number densities of galaxies

In this work, we have constrained the parameters σ₈, Ω_m and H₀ using galaxy samples with a number density of 3 × 10⁻³h³Mpc⁻³ at a redshift of z ∼ 0.1. We have verified that the theoretical framework successfully predicts both the Void Probability Function and the number density of voids larger than r for this number density. However, it remains to be checked whether the theoretical framework continues to successfully predict these two statistics at a different number densities of galaxies.

To check whether the theoretical framework works correctly for different number densities of galaxies we need to construct different samples (from Uchuu galaxy simulation box, for example) with different galaxy number densities, calculate the VPF and n_v(r) from that sample and compare with the predicted value given by the theoretical framework with the same coefficients that we have calculated in this work. The only parameter that must be changed for each sample is the mass, m_g. The number densities we use are n_g = {1 × 10⁻³, 2 × 10⁻³, 3 × 10⁻³, 4 × 10⁻³, 5 × 10⁻³, 6 × 10⁻³}h³Mpc⁻³.

In Figure H.1, we can observe the VPF (left panel) and number density of voids larger than r (right panel) predicted by theoretical framework with continuous lines and obtained in Uchuu galaxy simulation boxes for these galaxy number densities with dots and the abundance of voids larger than r (right panel). It can be seen that the theoretical framework predicts successfully the VPF for all galaxy number densities without any need to change the values of α or μ, except for galaxy number densities smaller than n_g ≤ 2 × 10⁻³h³Mpc⁻³ (see the irregularity obtained for r between 11 and 14 h⁻¹Mpc). Therefore, the theoretical framework is only valid for large galaxy number densities (n_g ≥ 2 × 10⁻³h³Mpc⁻³), and the coefficients α and μ depend on this galaxy number density.

Fig. H.1. VPF (left panel) and number density of voids larger than r (right panel) in real space predicted by the theoretical framework developed in this work (continuous lines) and measured in Uchuu galaxy box (points) for different galaxy number densities, ng (in units of h3Mpc−3), at redshift z = 0.092.

Appendix I: Content of the void catalogues used in this work

The columns of the void catalogues of halo and galaxy simulation boxes are the following

• X[MPC/H]: x-position of the centre of the void (comoving h⁻¹Mpc).

• Y[MPC/H]: y-position of the centre of the void (comoving h⁻¹Mpc).

• Z[MPC/H]: z-position of the centre of the void (comoving h⁻¹Mpc).

• RADIUS[MPC/H]: radius of the void (comoving h⁻¹Mpc).

For Uchuu-SDSS lightcones and SDSS, there are four additional columns:

• RA[DEG]: right ascension (degrees).

• DEC[DEG]: declination (degrees).

• ZOBS: observed redshift of the centre of the void (accounting for peculiar velocities. The fiducial cosmology used for SDSS voids is Planck 2015 (Planck Collaboration XIII 2016)).

• completeness: mean completeness of the void, calculated as the mean completeness of all the points uniformly distributed in its volume.

All Tables

All Figures

	Fig. 1. Halo mass function multiplied by the mean halo mass within virial radius (including unbound particles) of each bin of the four box catalogues used in this work. Shaded regions represent poissonian errors.
In the text

	Fig. 2. Voids with r > 9 h−1 Mpc (spheres) found in a region of P18 box catalogue (black and red points) with number density 3 × 10−3 h3 Mpc−3. Points that define the voids (i.e. those lying in their surface) are highlighted with a red circle. The volume of the box is 803 h−3 Mpc3.
In the text

Fig. 3. Number density of voids larger than r, nv(r)c, obtained for voids found in the distribution of SDSS (black points) and Uchuu-SDSS galaxies (pink solid line) with Mr < −20.5 and 0.02 < z < 0.132 and predicted by the theoretical framework described in Section 4 with Planck 2015 parameters (black dashed line). Shaded pink region delimits the standard deviation (1σ) of the 32 Uchuu-SDSS lightcones.

In the text

Fig. 4. Number density of voids larger than r in real space is shown for the theoretical framework (lines) and simulations (dots), shown in the left panel. The ratio between simulations and theoretical framework is shown in the bottom panels for the Uchuu, P18, Low, and VeryLow box catalogues with a number density of $\overline{n}=3\times {10}^{-3}{h}^3$ Mpc−3. The shaded region in the bottom panel delimits the region between 0.9 and 1.10 for the ratio.

In the text

Fig. 5. Number density of voids larger than r in redshift space is shown for the theoretical framework (lines) and simulations (dots) in the left panel, while the ratio between simulations and theoretical framework is shown in the bottom panels for the Uchuu, P18, Low, and VeryLow box catalogues with a number density of $\overline{n}=3\times {10}^{-3}{h}^3$ Mpc−3. The shaded region in the bottom panel indicates delimits the region between 0.9 and 1.10 for the ratio.

In the text

	Fig. 6. Marginalised posteriors in σ8 − Γ plane obtained for Uchuu-SDSS voids and SDSS voids voids using the maximum likelihood test with Bayesian approach with the theoretical framework developed in this work. We also show the marginalised posteriors for Planck 2018. The contours indicate the 68% (1σ) and 95% (2σ) credible intervals.
In the text

	Fig. 7. Marginalised posteriors in Ωm-H0, σ8-H0, and σ8 − Ωm planes obtained for SDSS redshift survey (left) and Uchuu-SDSS lightcones (right) using the maximum likelihood test with Bayesian approach with the theoretical framework developed in this work. The contours indicate the 68% (1σ) and 95% (2σ) credible intervals.
In the text

	Fig. 8. KiDS-1000, DESY3, and KiDS-1000+DESY3 contours, shown in the left, middle, and right panels, respectively. The SDSS voids contours are given in the σ8 − Ωm plane, using the same parameter ranges to make a fair comparison, along with the combination of SDSS voids with each weak lensing work. In Table 3, the parameter ranges of each work are indicated.
In the text

Fig. 9. Comparison between recent constraints on the parameters σ8 (left panel) and H0 (right panel) from different cosmological probes. The error bars represent 68% confidence intervals. The black error bars represent the values constrained in this work. The references of the rest of the works, from top to bottom in left panel are: Dark Energy Survey and Kilo-Degree Survey Collaboration (2023), Heymans et al. (2021), and Planck Collaboration VI (2020), and for the right panel: Freedman et al. (2020), Riess et al. (2022), Brout et al. (2022a,b), and Planck Collaboration VI (2020).

In the text

	Fig. 10. Comparison between recent constraints on S8 from different cosmological probes. The error bars represent 68% confidence intervals. The black error bars represent the values constrained in this work. The references of the rest of the works, from top to bottom are: Dark Energy Survey and Kilo-Degree Survey Collaboration (2023), Li et al. (2023), and Planck Collaboration VI (2020).
In the text

	Fig. A.1. Void identification process.
In the text

Fig. D.1. VPF in real space is shown for the theoretical framework (lines) and simulations (dots) given in the top panel, while the ratio between simulations and theoretical framework is shown in bottom panels for the Uchuu, P18, Low, and VeryLow box catalogues with number density $\overline{n}=3\times {10}^{-3}{h}^3$Mpc−3. Shaded region in bottom panel indicates delimits the region between 0.9 and 1.10 for the ratio.

In the text

Fig. D.2. VPF in redshift space is shown for the theoretical framework (lines) and simulations (dots) given in the top panel, while the ratio between simulations and theoretical framework is shown in bottom panels for the Uchuu, P18, Low and VeryLow box catalogues with number density $\overline{n}=3\times {10}^{-3}{h}^3$Mpc−3. Shaded region in bottom panel indicates delimits the region between 0.9 and 1.10 for the ratio.

In the text

	Fig. F.1. Constraints from Nv, i identified in Uchuu (first row), P18 (second row), Low (third row) and VeryLow (row) halo simulation boxes in real (first column) and redshift (second column) spaces using Maximum Likelihood test. The contours indicate the 68% (1σ) and 95% (2σ) credible intervals. The black dots are the real values of σ8 and H0 of each simulation.
In the text

	Fig. G.1. VPF (left panel) and number density of voids larger than r (right panel) predicted by the theoretical framework developed in this work (continuous lines) and measured in Uchuu galaxy box (points) for different galaxy number densities, ng, at redshift z = 0.092.
In the text

	Fig. H.1. VPF (left panel) and number density of voids larger than r (right panel) in real space predicted by the theoretical framework developed in this work (continuous lines) and measured in Uchuu galaxy box (points) for different galaxy number densities, ng (in units of h3Mpc−3), at redshift z = 0.092.
In the text