The abundance of voids and the excursion set formalism

Abstract

We present measurements of the number density of voids in the dark matter distribution from a series of N-body simulations of a Λ cold dark matter cosmology. We define voids as spherical regions of ρ_v = 0.2ρ_m around density minima in order to relate our results to the predicted abundances using the excursion set formalism. Using a linear underdensity of δ_v = −2.7, from a spherical evolution model, we find that a volume-conserving model, which does not conserve number density in the mapping from the linear to non-linear regime, matches the measured abundance to within 16 per cent for a range of void radii 1 < r(h⁻¹ Mpc) < 15. This model fixes the volume fraction of the universe which is in voids and assumes that voids of a similar size merge as they expand by a factor of 1.7 to achieve a non-linear density of ρ_v = 0.2ρ_m today. We find that the model of Sheth and van de Weygaert for the number density of voids greatly overpredicts the abundances over the same range of scales. We find that the volume-conserving model works well at matching the number density of voids measured from the simulations at higher redshifts, z = 0.5 and 1, as well as correctly predicting the abundances to within 25 per cent in a simulation of a matter dominated Ω_m = 1 universe. We examine the abundance of voids in the halo distribution and find fewer small, r < 10 h⁻¹ Mpc, voids and many more large, r > 10 h⁻¹ Mpc, voids compared to the dark matter. These results indicate that voids identified in the halo or galaxy distribution are related to the underlying void distribution in the dark matter in a complicated way which merits further study if voids are to be used as a precision probe of cosmology.

Cosmology: theory, large-scale structure of Universe

1 INTRODUCTION

Galaxy redshift surveys allow us to study and map out the large-scale structure of our Universe revealing a hierarchical mass distribution with substructure over a wide range of scales. The main components of the galaxy distribution are arranged in a remarkable ‘cosmic web’ (Bond, Kofman & Pogosyan 1996) made up of clusters of galaxies connected by filaments with large empty voids which occupy most of the volume. Only recently have systematic studies using voids as precision probes of the growth of structure been possible due to the increased depth and volume of current galaxy surveys (York & SDSS Collaboration 2000; Colless et al. 2001; Abazajian et al. 2009). In this paper, we study the distribution of underdense void regions in the dark matter and halo distributions using N-body simulations. We focus on the excursion set method which gives an analytical prescription for the number density of voids which we compare with measurements from simulations.

Voids are a common feature in galaxy surveys with one of the most well-known discoveries being the void in Boötes which has a diameter of ∼50 h⁻¹ Mpc (Kirshner et al. 1981). Since then several surveys such as the Center for Astrophysics Redshift Survey (Geller & Huchra 1989), the Southern Sky Redshift Survey (Maurogordato, Schaeffer & da Costa 1992) and the deeper Las Campanas Redshift Survey (Shectman et al. 1996) have identified voids in the distribution of galaxies and clusters confirming that they are the dominant and volume filling component of our Universe. Most recently, Pan et al. (2012) and Sutter et al. (2012a) both used the Sloan Digital Sky Survey Data Release 7 (SDSS DR7) (Abazajian et al. 2009) to identify voids. Pan et al. (2012) found 1054 statistically significant voids with radii r > 10 h⁻¹ Mpc with an absolute magnitude cut of M_r < −20.09. They argue that voids of effective radius r_eff ∼ 20 h⁻¹ Mpc dominate the void volume with the largest void in the sample having r ∼ 30 h⁻¹ Mpc. Sutter et al. (2012a) constructed the first public void catalogue using the full extent of the SDSS DR7 spectroscopic survey which included the LRGs and found large voids in the sample of r ∼ 50–60 h⁻¹ Mpc in radius.

Early numerical and theoretical work on the evolution of voids by Regos & Geller (1991), Blumenthal et al. (1992) and Dubinski et al. (1993) focused on the expansion of initial linear undensities up to the moment of shell crossing, which is used to define a characteristic time in the formation of voids. They considered spherical voids in a Ω_m = 1 universe and found that shell crossing occurs at a linear underdensity of −2.7 at which point the comoving size of the void has increased by a factor of 1.7 (see also van de Weygaert & van Kampen 1993; Friedmann & Piran 2001). Many studies since then have focused on analysing the dynamics and statistical properties of voids such as the void probability function (VPF), the probability that a randomly placed sphere will contain no objects (White 1979), the void filling factor, the void number density and void density profiles (see e.g. van de Weygaert 1991; Mathis & White 2002; Benson et al. 2003; Colberg et al. 2005; Shandarin et al. 2006; Betancort-Rijo et al. 2009; Einasto et al. 2011; van de Weygaert & Platen 2011; Kreckel et al. 2012; Aragon-Calvo & Szalay 2013).

Recent studies have looked at stacking voids in order to increase the statistical significance of weak lensing signals (Higuchi, Oguri & Hamana 2012; Krause et al. 2013) (see also Amendola, Frieman & Waga 1999), the Integrated Sach–Wolfe effect (Cai et al. 2013; Ilic, Langer & Douspis 2013) or to extract cosmological parameters by modelling the distortions in redshift space (Bos et al. 2012; Lavaux & Wandelt 2012; Sutter et al. 2012b) or as a test of modified gravity (Clampitt, Cai & Li 2013). The precision of these tests relies on many factors, for example, given a survey or numerical simulation of a certain size, how robustly can we measure statistics for a void of a given size, how accurately can we predict the number density of voids in the galaxy/dark matter halo distribution and how well do voids in the galaxy/halo distribution trace voids in the dark matter. In this paper, we address these three issues. In our discussion of a robust void finder, we do not compare with all the other algorithms which have been used in previous studies – our choice of void finder is motivated by the excursion set formalism for the abundance of voids which we aim to test.

In analysing both galaxy surveys and numerical simulations a wide variety of void finding algorithms have been used to define underdense regions as voids. Colberg et al. (2008) carried out the first systematic review of 13 different void finders, identifying only two areas of agreement amongst the different algorithms: that voids are very underdense (ρ ∼ 0.05ρ_m) at their centres and that voids have very steep edges. The void finding methods include the construction of proto-voids around local minima in the smoothed density field, after separating the galaxy sample into ‘wall’ and ‘void’ galaxies (see e.g. El-Ad & Piran 1997; Hoyle & Vogeley 2002); merging proto-voids which results in non-spherical voids (Colberg et al. 2005); the watershed algorithm (Platen, van de Weygaert & Jones 2007) which uses the dtfe method (Schaap 2007; Cautun & van de Weygaert 2011). The watershed void finder identifies minima in the density field and constructs voids by flooding basins until the ‘landscape’ resembles a segmented plane where the edges of each segment outline a void region. A similar algorithm, which we make use of in this paper, is the zobov (Neyrinck 2008) void finder which uses the Voronoi tessellation field method (see e.g. van de Weygaert 2007) to partition particles into zones, which are then joined together near density minima, into voids.

There have been many studies of the excursion set method to predict the abundance of dark matter haloes in the Universe (see e.g. Zentner 2007, for a review). In comparison, fewer studies have focused on testing the excursion set predictions for underdensities in the dark matter distribution and we briefly outline some of these works here. In applying this method to voids, Sheth & van de Weygaert (2004) presented a model for the abundance of voids in the dark matter which included the influence of the larger scale environment on the formation of a void. Their model takes into account two effects, first, a void of a given size may be embedded in another underdense region which is on a larger scale, the ‘void-in-void’ scenario, and secondly, a void of a given size could be embedded in an overdense region on a larger scale, the ‘void-in-cloud’ scenario. Furlanetto & Piran (2006) analysed how the barrier for shell crossing of a void in the galaxy distribution would differ from the linear theory barrier for dark matter, finding that voids selected from catalogues of luminous galaxies should be larger than those selected from faint galaxies (see also D'Aloisio & Furlanetto 2007). D'Amico et al. (2011) consider using voids as a probe of primordial non-Gaussianity and calculate the abundance of voids using the excursion set formalism and the two barrier prescription of Sheth & van de Weygaert (2004). Shandarin et al. (2006) define voids as isolated regions of the low-density excursion set specified by density thresholds and measured the abundance and morphology of voids using N-body simulations. In this paper, we wish to test the predictions of the Sheth & van de Weygaert (2004) excursion set model by comparing them to measurements of void abundances from N-body simulations. Our definition of a void is similar to that adopted by Shandarin et al. (2006) as we use a strict density threshold to define the void edge (although we do not use isodensity contours). To our knowledge, this is the first time that this model has been directly compared with numerical simulations.

This paper is organized as follows. In Section 2, we discuss the excursion set method as it applies to dark matter haloes and voids. Appendices A and B review the salient features of the spherical evolution model that connects the two. In Section 3, we detail our void finder and the N-body simulations that were carried out. In Section 4, we present the main results of this paper on the number density of voids in three different cosmological models at z = 0 and we show how this abundance changes with redshift. We also present the measured number density of voids in the halo distribution. In Section 5, we present our conclusions.

2 EXCURSION SET FORMALISM FOR VOID ABUNDANCE

In this section, we begin with a brief review of the excursion set formalism in Section 2.1. It is well known that in combination with spherical collapse this approach provides insight into many aspects of halo formation and can be used to predict dark matter halo abundances and clustering (see e.g Zentner 2007, for a review of the subject). The analogous spherical expansion model can likewise be used to make excursion set predictions for voids (Sheth & van de Weygaert 2004). We review this extension in Section 2.2 and show that it requires modifications on physical grounds. We propose a simple modification based on volume fraction conservation in Section 2.3.

2.1 Excursion set formalism

The excursion set formalism at its heart relies on knowledge of the statistical properties of the linear density field. In Fourier space, the linear density fluctuation field smoothed on a scale R is given by

\begin{eqnarray} \delta (\boldsymbol {x},R) = \int \frac{{\rm d}^3k}{(2\pi )^3} \delta (\boldsymbol {k}) W(\boldsymbol {k},R) {\rm e}^{-{\rm i}\boldsymbol {k} \cdot \boldsymbol {x}} \,, \end{eqnarray}

(1)

where |$\delta (\boldsymbol {k})$| is the Fourier transform of the density perturbation |$\delta (\boldsymbol {x}) = [\rho (\boldsymbol {x}) - \rho _{\rm m}]/\rho _{\rm m}$|⁠, |$\rho (\boldsymbol {x})$| is the local density at comoving position |$\boldsymbol {x}$|⁠, ρ_m is the background matter density and |$W(\boldsymbol {k},R)$| is a filter function in Fourier space. It is common to relate the smoothing scale R to the corresponding variance of the linear density field

\begin{eqnarray} \sigma ^2(R) \equiv S(R)= \int \frac{{\rm d}k}{k} \, \frac{k^3 P(k)}{2\pi ^2} |W(k,R)|^2 \,, \end{eqnarray}

(2)

where P(k) is the matter power spectrum in linear perturbation theory. We can refer to a trajectory |$\delta (\boldsymbol {x},S)$| as a sequence of overdensities given by subsequent increases in the smoothing scale by increments ΔS. When a tophat filter in k-space is used then |$\delta (\boldsymbol {x},S)$| executes a random walk. Given an underlying Gaussian distribution for the linear density field, the excursion set formalism allows us to associate probabilities to random walks that satisfy a given set of criteria for the smoothing scale at which they cross various density thresholds. Its use in defining the statistics of objects in the non-linear regime requires a model that associates such criteria to objects.

2.2 Spherical evolution and SVdW model

The spherical evolution model provides a complete description of the non-linear evolution of a spherically symmetric top-hat density perturbation. One of the main features of this model is that the evolution does not depend on the initial size of the region, i.e. on the initial radius or enclosed mass, but only on the amplitude of the initial top-hat overdensity.

For the collapse of perturbations, the spherical evolution model in combination with the excursion set provides a good description of the statistics of dark matter haloes. As we review in Appendix A, collapse occurs when the linear density fluctuation reaches a critical value or barrier δ_c. We can then use the excursion set formalism to determine the fraction of trajectories that cross this barrier for the first time, accounting for the cloud-in-cloud process, within some dln σ of a smoothing scale σ through the differential fraction

\begin{eqnarray} f_{\ln \sigma }(\sigma ) \equiv \frac{{\rm d} f}{{\rm d }\ln \sigma }= \sqrt{\frac{2}{\pi }} \frac{\delta _{\rm c}}{\sigma } {\rm e}^{-\frac{\delta ^2_{\rm c}}{2\sigma ^2}} \,. \end{eqnarray}

(3)

Since both mass and number are conserved in the collapse, the linear theory mapping σ(M) carries over to the non-linear regime and so the mass function, or the comoving differential number density of haloes is

\begin{eqnarray} \frac{{\rm d}n}{{\rm d}\ln M} = \frac{\rho _{\rm m}}{M} f_{\ln \sigma }(\sigma ) \frac{{\rm d }\ln \sigma ^{-1}}{{\rm d }\ln M}, \end{eqnarray}

(4)

where ρ_m/M is the number density of such objects if the fraction were unity.

We can extend the model to underdense regions in the initial density field. These are naturally associated with voids in the evolved density field today. A key assumption in making the connection between the excursion set and the abundance of non-linear objects is that each collapse occurs in isolation. This makes sense for collapsing objects since the comoving volume occupied shrinks. In contrast to overdense regions which contract, voids expand. We shall see that this causes a problem for mapping excursion set predictions on to the statistics of voids.

Nonetheless, let us start with the simple spherical evolution model following Sheth & van de Weygaert (2004). The critical density threshold is defined to be when the expanding shells cross (see e.g. Fillmore & Goldreich 1984; Suto, Sato & Sato 1984; Bertschinger 1985). As shown in Appendix A for an Einstein–de Sitter (EdS) universe, this occurs when the non-linear average density within the void reaches ρ_v = 0.2ρ_m or when the linear density threshold reaches δ_v = −2.7. Note that we will use this notation of ρ_v to refer to the non-linear density of the void region and δ_v to refer to the linear underdensity used as a threshold in the excursion set model. We show in Appendix A that these EdS values suffice for the accuracy to which we wish to describe alternate cosmologies such as the Λ cold dark matter (ΛCDM) model.

Once we have this value for the void barrier, we can follow the excursion set formalism for determining the fraction of random walks which pierce the barrier δ_v. Similar to the cloud-in-cloud process, the void-in-void process accounts for the fact that a void of a given size may be embedded in another underdense region on a larger scale. We thus define the first crossing distribution by associating the random walks with the smoothing scale for which they first cross the barrier δ_v.

The second process, the void-in-cloud scenario, occurs when a void of a given size is embedded in an overdense region on a larger scale, which will eventually collapse to a halo and squash the void out of existence. In order to account for the void-in-cloud effect, Sheth & van de Weygaert (2004) proposed that the excursion set method applied to voids requires a second barrier, the threshold for collapse of overdense regions, δ_c. In calculating the first crossing distribution, Sheth & van de Weygaert (2004) argued that we need to determine the largest scale at which a trajectory crosses the barrier δ_v given that it has not crossed δ_c on any larger scale. They posit that the value of δ_c should lie somewhere in between δ_c = 1.06, the value at turnaround in the spherical collapse model, and δ_c = 1.686, the value at the point of collapse (see also Paranjape, Lam & Sheth 2012). In this paper, we shall refer to the model of Sheth & van de Weygaert (2004) as the ‘SVdW’ model.

By the same reasoning as applied to haloes, the SVdW formula for the abundance of voids in linear theory is given by

\begin{eqnarray} \frac{{\rm d} n_{\rm L} }{{\rm d }\ln M} = \frac{\rho _{\rm m}}{M} f_{\ln \sigma }(\sigma ) \frac{{\rm d }\ln \sigma ^{-1}}{{\rm d }\ln M} \,, \end{eqnarray}

(5)

where

\begin{eqnarray} f_{\ln \sigma }(\sigma ) &=&2 \sum _{j=1}^\infty {\rm e}^{-\frac{(j \pi x)^2}{2}} j\pi x^2 \sin (j\pi \mathcal {D}), \end{eqnarray}

(6)

with

\begin{eqnarray} \mathcal {D} = \frac{|\delta _{\rm v}|}{\delta _{\rm c}+ |\delta _{\rm v}|} \,,\quad x =\frac{ \mathcal {D} }{|\delta _{\rm v}|}\sigma \,. \end{eqnarray}

(7)

Note that Sheth & van de Weygaert (2004) give f_ln S = Sdf/dS = f_ln σ/2. We have added the subscript ‘L to remind the reader that the logic relies on equating a number density derived from linear theory to the number density of some non-linear object for reasons that will be clear below.

Since the infinite series in equation (6) is cumbersome to work with, it is useful to have an accurate closed form expression. As we discuss in Appendix B, the accuracy of the approximation given in Sheth & van de Weygaert (2004) is uncontrolled as σ → ∞. Instead, we find the limiting forms for equation (6) such that the domain of validity of the approximation is well defined. Note that as σ → 0, the opposing barriers are high and the sum must return the single barrier expression since the probability of first crossing the collapse barrier is vanishingly small. This fixes the form as x → 0. As |$\mathcal {D}$| increases towards unity, we lower the collapse barrier relative to the void barrier and the value of x at which this limit is approached decreases. Correspondingly to achieve a matching at this point, we need to keep more terms in the sum. The largest value that we will be interested in is |$\mathcal {D}<3/4$| and so it suffices to keep four terms

\begin{eqnarray} f_{\ln \sigma }(\sigma ) &\approx & \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\sqrt{\frac{2}{\pi }}\frac {|\delta _{\rm v}|}{\sigma } {\rm e}^{-\frac{\delta _{\rm v}^2}{2\sigma ^2} },&x\le 0.276 \\ 2 \sum _{j=1}^4 {\rm e}^{-\frac{(j \pi x)^2}{2}} j\pi x^2 \sin (j\pi \mathcal {D}), & x>0.276\\ \end{array}\right. \end{eqnarray}

(8)

which is accurate at the 0.2 per cent level or better across the domain of validity. The approximation of equation (8) is used in all the numerical work throughout the paper.

We can alternately express the number density in terms of the linear theory radius of the void r_L. Using ρ_m/M = 1/V(r_L) and defining the volume of a spherical region of an arbitrary radius, R, as

\begin{eqnarray} V(R) \equiv \frac{4}{3} \pi R^3 \,, \end{eqnarray}

(9)

we obtain

\begin{eqnarray} \frac{{\rm d} n_{\rm L} }{{\rm d }\ln r_{\rm L}} = \frac{ f_{\ln \sigma }(\sigma ) }{V(r_{\rm L})} \frac{{\rm d }\ln \sigma ^{-1}}{{\rm d }\ln r_{\rm L}} \,. \end{eqnarray}

(10)

In the spherical evolution model, the actual void expands from its linear radius. At the epoch of shell crossing ρ_v = 0.2ρ_m. Given that

\begin{eqnarray} \frac{r}{r_{\rm L}} = \left( \frac{\rho _{\rm m}}{\rho _{\rm v}} \right)^{1/3} \,, \end{eqnarray}

(11)

spherical expansion predicts that this expansion factor is r ≈ 1.7r_L. The void abundance therefore becomes

\begin{eqnarray} \frac{{\rm d} n }{{\rm d }\ln r} = \frac{{\rm d} n_{\rm L} }{{\rm d }\ln r_{\rm L}} \Big |_{r_{\rm L}= r/1.7}\,. \end{eqnarray}

(12)

Note that in this model dn/dln r shifts left to right in scale through the non-linear growth but does not change in amplitude, as is shown in Fig. 1.

Figure 1.

Void abundance model predictions. In the SVdW model, the number density of linear underdensities (blue curve) remains unchanged in void formation and only their sizes change (arrow to orange curve). In the V dn model, the number density also changes so as to conserve the volume fraction in voids, lowering the amplitude at fixed shape (arrow to grey curve). Varying 1.06 ≤ δ_c ≤ 1.686 (shaded or hatched regions) changes the abundance significantly only for small voids r ≲ 1 h^{− 1} Mpc. We take δ_v = −2.7 throughout. We use the σ₈ = 0.8 ΛCDM cosmology as listed in Table 2 here and in the following figures unless otherwise stated.

The SVdW model has two parameters δ_c and δ_v. The latter is fixed by the shell-crossing criterion whereas the former is expected to vary within 1.06 ≤ δ_c ≤ 1.686. In Fig. 1, we also show that for the range of radius of interest (r > 1 h^{− 1}Mpc), changing δ_c within its expected range has little effect on the void abundance.

The SVdW model makes a very specific prediction for the abundance of large voids. Again the key assumption of the SVdW model is that the comoving number density of objects is conserved during the evolution n = n_L and only their size has changed. Unfortunately, for spherical evolution this assumption is invalid for large voids. In particular, the cumulative volume fraction in voids larger than R defined as

\begin{eqnarray} \mathcal {F}(R) = \int _R^\infty \frac{{\rm d }r}{r} \, V(r) \frac{{\rm d}n}{{\rm d}\ln r}, \end{eqnarray}

(13)

exceeds unity for radii of interest. In Fig. 2, we demonstrate that this problem cannot be cured by changing δ_c within the expected range as it only affects small voids whereas the problem appears at R ≈ 2 h^{− 1}Mpc. Indeed, if we take R → 0, then for the exact f_ln σ given by equation (6) (Sheth & van de Weygaert 2004)

\begin{eqnarray} \mathcal {F}(0) = \left( \frac{r}{r_{\rm L}} \right)^3 \int _0^\infty \frac{d \sigma }{\sigma } f_{\ln \sigma } = \left( \frac{r}{r_{\rm L}} \right)^3 (1-\mathcal {D}). \end{eqnarray}

(14)

This result suggests that reducing δ_v → 0 simultaneously takes r → r_L and |$\mathcal {D} \rightarrow 0$| bringing SVdW asymptotically to physicality |$\mathcal {F}(0)=1$|⁠. Strictly speaking, δ_v is fixed by the shell-crossing criterion. However, given the approximate nature of the correspondence between the isolated spherical expansion model and real voids, it is interesting to explore whether modifications to this criterion can bring the SVdW model into agreement with physicality and simulations. If we change the non-linear density at which voids are defined ρ_v/ρ_m, the linear density threshold δ_v and the expansion factor r/r_L must change in a self-consistent fashion (see Fig. A1 and equation B1). In Fig. 3, we show that changing δ_v alters the shape of the abundance function. As |δ_v| decreases, the steepness of the abundance function also decreases. Thus, although lowering δ_v can make the total volume fraction physical (Fig. 3, lower panel), it increases the abundance of the largest voids. We shall see that the agreement between simulations and the abundance of voids in the excursion set method is remarkably good if we remove the assumption of isolated spherical expansion. Assuming that the number density of voids is conserved as they expand causes the SVdW model to greatly overpredict the abundance of large voids regardless of the choice of δ_c and δ_v.

Figure 2.

The cumulative volume fraction in voids with radii larger than R for the various models: linear theory (blue striped region, R = r_L), SVdW model (orange striped region, R = r), V dn model (grey shaded region, R = r). Regions correspond to the expected range of 1.06 ≤ δ_c ≤ 1.686 and we take δ_v = −2.7 throughout. For SVdW, the fraction unphysically exceeds unity at R ≈ 2 h^{− 1}Mpc while for the V dn model conserves the total fraction from the linear theory of |$\mathcal {F}(0)\approx 0.3$|⁠.

Figure 3.

Relaxing the shell-crossing criterion of the SVdW model void abundance predictions. Upper: variation of the void underdensity ρ_v/ρ_m changes both the shape of the abundance through the linear barrier δ_v and the size of the voids or horizontal shift through r/r_L = (ρ_v/ρ_m)^−1/3. Decreasing |δ_v| increases the number of large voids and decreases that of small voids. Lower: the cumulative volume fraction in voids with radii larger than R decreases as |δ_v| → 0 and R → 0 but at the expense of making the larger voids more abundant. We use δ_c = 1.06 throughout.

2.3 Volume-conserving V dn model

We propose a simple fix to the unphysicality of the isolated spherical expansion model for voids. We require that the volume fraction and shape of the abundance function is fixed during the expansion, rather than assuming that the expansion of isolated voids preserves their total number density. Specifically, if we define the volume fraction in linear theory, |$\mathcal {F}_{\rm L}$|⁠, as

\begin{eqnarray} \mathcal {F}_{\rm L}(R_{\rm L}) = \int _{R_{\rm L}}^\infty \frac{ {\rm d } r_{\rm L}}{r_{\rm L}} V(r_{\rm L}) \frac{{\rm d}n_{\rm L}}{{\rm d} \ln r_{\rm L}}\,, \end{eqnarray}

(15)

then this fraction is conserved if we define the non-linear abundance as

\begin{eqnarray} V(r)\, {\rm d}n = V(r_{\rm L}) {\rm d}n_{\rm L}\big |_{r_{\rm L}(r)}. \end{eqnarray}

(16)

In this picture, when a void expands from r_L → r it combines with its neighbours to conserve volume and not number. Thus, the abundance becomes

\begin{eqnarray} \frac{{\rm d} n }{{\rm d }\ln r} &=& \frac{ V(r_{\rm L}) }{V(r)} \frac{{\rm d} n_{\rm L}}{{\rm d }\ln r_{\rm L}} \frac{{\rm d} \ln r_{\rm L}}{{\rm d }\ln r} \Big |_{r_{\rm L}(r)} \nonumber \\ &=& \frac{ f_{\ln \sigma }(\sigma )}{V(r)} \frac{{\rm d }\ln \sigma ^{-1}}{{\rm d }\ln r_{\rm L}} \frac{{\rm d} \ln r_{\rm L}}{{\rm d }\ln r} \Big |_{r_{\rm L}(r)} \,. \end{eqnarray}

(17)

We call this model the V dn model and show its abundance prediction in Fig. 1. We have left the mapping r(r_L) general here since the specific form from isolated spherical expansion until shell crossing may not apply here. We will however adopt r = 1.7r_L for voids with non-linear density ρ_v = 0.2ρ_m from N-body simulations as a starting point. Note that in this case dln r_L/dln r = 1 and the impact of going from the linear to the non-linear abundance is both a shift in scale and a change in amplitude with no change in shape, as is shown as the combination of arrows in Fig. 1.

In Fig. 2, we also show the cumulative volume fraction with this abundance function, along with that for linear theory defined in equation (15). Since by construction the volume fraction is conserved, the two curves differ only by a horizontal shift in scale.

Since the V dn model is not the unique means of constructing a physical model, it is interesting to explore other ways of keeping the volume fraction below unity. Phenomenologically, we can decouple the relationship in equation (B1) between the parameters δ_v and r/r_L provided by the spherical expansion model. In fact, we can choose these parameters so as to mimic the V dn predictions for a fixed cosmology. For example, in the upper panel of Fig. 4, we find we can change the parameters δ_v → −2 and r/r_L → 1 in the SVdW model to fit the V dn model in the σ₈ = 0.8 ΛCDM cosmology listed in Table 2. However, this change then predicts very different abundances than the V dn model for a different cosmology as shown with the EdS cosmology listed in Table 2 and Fig. 4 (lower panel). We shall show below that simulation results favour the V dn model over universal changes in δ_v and r/r_L. The V dn model retains the δ_c parameter from the SVdW model to describe the influence of surrounding mass concentrations on the growth of voids. Within the excursion set method this parameter accounts for the crushing of small voids which reside in overdense regions. Note that the assumption of spherical expansion of these small voids, as well as the spherical collapse of the larger overdense region should be taken as a simple approximation which will not be accurate for small non-spherical voids. In this study, we focus on testing the model using simulations of voids with radii r > 1 h⁻¹ Mpc whose abundance is not affected by this crushing effect.

Figure 4.

Void abundance in the V dn model (black dot–dashed curve) and a modified SVdW model (blue solid line) with ad hoc variations designed to fit the ΛCDM V dn model. Upper panel: we choose δ_v = −2 and r/r_L = 1 in the SVdW model in violation of spherical expansion predictions in a σ₈ = 0.8 ΛCDM cosmology. Lower panel: we show that the same set of parameters give a poor fit in an EdS model (see Table 2). For all curves we use δ_c = 1.686.

3 SIMPLE VOID FINDING ALGORITHM

In Section 3.1, we outline the void finding algorithm used to identify voids in both the dark matter and halo populations in this work. In Section 3.2, we present the details of the N-body simulations which were carried out as summarized in Table 2.

3.1 Void finder

As we have already mentioned one of the main complications in studying the distribution of voids in the large-scale structure of the Universe is finding a robust definition of what a void is (see Colberg et al. 2008, for a comparison of void finders). In this work, we wish to make a direct comparison to the predictions of the excursion set formalism which assumes that these underdense regions are non-overlapping spheres of a given underdensity corresponding to a region at the moment of shell crossing. We shall retain the moment of shell crossing as the key feature which defines a non-linear void in the matter distribution today although we also compare the measured abundance of voids with different underdensities to the predictions of the volume-conserving model in Section 4.1.

We start with the publicly available code zobov (Neyrinck 2008) which uses Voronoi tessellation to estimate densities and find both voids and subvoids. The main advantage of using tessellation methods is that it gives a local density estimate by dividing space into cells, where the cell around any given particle is the region of space closer to that particle than to any other. The Voronoi tessellation also gives a natural set of neighbours for each particle which zobov uses to construct zones around density minima.

The output from zobov is useful for our purposes for two main reasons. First, it outputs a linked list of zones in the dark matter distribution, which is also ordered by density contrast. It thus provides a tree structure which we can prune according to the definition of a void. Secondly, zobov identifies the ‘core’ or least dense particle in a zone and returns its density as well as a measure of the probabilities that each collection of zones arises from Poisson fluctuations. Note that the list which zobov returns contains zones of various densities and Poisson probabilities, some of which could be overdense or not statistically significant, so it is necessary to prune the output from zobov in order to construct a void catalogue.

In constructing our void finder the goal is to identify all spherical non-overlapping underdense regions of average density ρ_v = 0.2ρ_m in a dark matter simulation. We use the output from zobov and find spherical regions centred around the core particle (lowest density particle) in a zone, which can encompass any particles which are around the zone returned by zobov and not necessarily part of the particular zone or collection of zones returned by zobov.

One of the outputs from zobov is a text file which lists individual zones and joined zones which are added to the list in a sequential process analogous to water flooding a plane with troughs of various heights (see also Platen et al. 2007). During flooding, when water from a particular zone or joined zones flows into a neighbouring deeper zone, the process stops and the zone is recorded in the list. A toy example of the output is shown in Table 1 where the zones are listed in order of density contrast. Fig. 5 shows an illustration of the Voronoi tessellation of the region surrounding Void #1 output from zobov and given in Table 1 which is made up of zones a (blue), b (red) and c (cyan). In this figure, the core particle (CorePar = 26) of zone a is shown as a green cross. Here, FileVoid# and CorePar refer to unique identification tags for the void and its core particle, respectively. CoreDen is the density, in units of the mean, of the void's core particle.

Figure 5.

An illustration of a spherical void identified using the zones output from zobov after Voronoi tessellation of the region. The Void #1 output from zobov given in Table 1 is shown as three shaded zones, a (blue), b (red) and c (cyan). The core particle of zone a is shown as a green cross while the void we identify in this region of a given density is shown as a green circle.

Table 1.

A toy example showing the first five columns from a zobov output file.

Void#(zones)	FileVoid#	CorePar	CoreDen	ZoneVol
1 (a, b, c)	2945	26	1.08e−02	4.9e+03
2 (b, c)	5033	83	1.8e−02	7.8e+02
3 (c)	1814	45	2.0e−01	1.9e+02

In order to count non-overlapping regions with an average density of 0.2, we perform the following two stages of analysis on the output from zobov. First, starting from the top of the zobov output file we determine if the first collection of zones listed e.g Void #1, which is made up of zones a, b and c, pass the following criteria:

• [r_min, r_max]: the radius corresponding to a sphere of equal volume should be > r_min and < r_max.

• The core particle density is <0.2ρ_m.

As we are searching for regions which have an average density ρ_v = 0.2ρ_m, we also only consider zones in the list which have a density ρ_v ≥ 0.2ρ_m in order to speed up the search.

If the collection of zones fulfils all of the above, then we proceed to the second stage. A spherical region around the centre is found by iteratively including one particle at a time moving away from the centre of the void until ρ_v = 0.2ρ_m. We assume that the volume corresponds to a sphere with radius equal to the distance from the centre to the last particle included. If for example, Void #1 did not pass the r_max, r_min criterion, then we consider if the deepest zone, zone a, does and if so we grow a sphere around the centre of this single zone. This step is important as the output from zobov only lists the deepest zone ‘a’ once and if Void #1 fails the r_max, r_min criteria the void finder would miss counting the deepest zone in the simulation box. In Fig. 5, we show an illustration of the spherical void (green circle) which is grown around the core particle (green cross) for Void #1. The spherical region is not necessarily restricted to include particles in the zones returned by zobov.

We then proceed to the next line in the output file and perform the same two stages of analysis. All of the particles in the spherical regions which are grown have been tagged and at any stage if there is any overlap of spheres, we disregard the less underdense zone to avoid double counting any volume in the simulation. In Fig. 5, this corresponds to growing another sphere around the core particle in zone b, until the required density is reached. This spherical void is added to the catalogue if it does not overlap with the void around zone a (green circle). Note that if we consider the output from zobov as a tree structure, then this procedure is similar to walking the tree from root to tip, pruning any branches after our criteria are met.

Using a cut in r_max and r_min as above allows us to avoid considering spuriously small voids and the first output in the text file which is a void which takes up the entire simulation box; however, we have checked that our results are not sensitive to changes in r_max and r_min but we retain these criteria in the void finder to speed up computation. For our simulations, we use the following: [r_min, r_max] = [0.5, 15], [1, 30], [2, 60] and [4,120] h⁻¹ Mpc for the 64, 128 and 256 and 500 h⁻¹ Mpc boxes, respectively.

We tested several different criteria in identifying voids and found that the two points listed above are sufficient to identify significant non-overlapping regions of a given underdensity in the particle distribution. In testing the robustness of our void finder, we considered the impact of the following adjustments to the method:

• As an alternative to using the core particle to define the centre of the void, we can use a volume-weighted centre defined as

  \begin{eqnarray} \boldsymbol {c} = \frac{\sum _i \boldsymbol {x}_iV_i}{\sum _i V_i}, \end{eqnarray}

  (18)

  where |$\boldsymbol {x}_i$| and V_i are the position and volume of each particle in the zones returned by zobov, respectively. We found that the abundance of voids was not substantially affected by this choice and so we use the core particle as the centre of the spherical region. Note that using the volume-weighted centre is more robust when using stacked voids (see e.g. Lavaux & Wandelt 2012).

• Instead of allowing the spherical region to include all particles which are around the void, we consider restricting it so that only particles which zobov list as being part of the zone are included when growing the sphere. We find that in the majority of cases the region of underdensity ρ_v = 0.2ρ_m is contained within the collection of zones returned by zobov and so this restriction does not affect the measured abundance of voids. Our results in Section 4 use all particles within 1.5 times the radius of the zones to find the spherical void. Including particles only within this radius was found to be sufficient considering the original collection of zones was required to have an average density of >0.2ρ_m. An alternative approach to this would be to use the actual volume of each zone particle when trying to find a void of a given average density. This would allow for irregularly shaped voids which it could be argued is a more ‘natural’ description of an actual void; however, as we mentioned we are trying to compare with the excursion set theory for abundances which assumes spherical voids.

• We originally included a third criterion in our void finder by requiring that the probability of a zone arising from a Poisson process was less than a given significance (see Neyrinck 2008, for more details). However, in the context of our void finder, we found that the core particle density requirement by itself was sufficient to get rid of spurious voids. This also agrees with the findings of Neyrinck (2008).

• In the algorithm we have described, we stop growing a sphere around the core particle when we find the desired underdensity at the maximum radius at which this occurs within the radius of the collection of zones. This is a different approach to simply stopping to record the first radius where ρ_v = 0.2ρ_m which would not take into account void-in-void scenarios. In practice, we find that accounting for a void-in-void effect alters the measured abundances by a small amount (e.g ∼7 per cent over the range 1 < r(h⁻¹ Mpc) < 10).

• The above method does not allow any overlap of voids within the simulation in order to compare with the excursion set method. In practice, for regions of ρ_v = 0.2ρ_m we found that the overlap was very small for the simulations we consider.

• zobov is run using all the particles in the simulation with a run-time density threshold parameter which can limit the growth of a collection of zones into high-density regions. We set this parameter to 0.2; however, we have verified that changing this run-time parameter has little effect on the abundance of underdensities found by our void finder.

3.2 N-body simulations

We measure the abundance of voids in the dark matter distribution using a series of N-body simulations in various box sizes. These simulations were carried out at the University of Chicago using the TreePM simulation code gadget-2 (Springel 2005). The ΛCDM model used has the following cosmological parameters: Ω_m = 0.26, Ω_DE = 0.74, Ω_b = 0.044, h = 0.715 and a spectral tilt of n_s = 0.96 (Sánchez et al. 2009). The linear theory rms fluctuation in spheres of radius 8 h⁻¹ Mpc is set to be σ₈ = 0.8 for our main simulation set of eight independent realizations of the ΛCDM cosmology. In order to investigate the abundance of voids in different cosmologies, we also carry out two additional simulations; one with a ΛCDM cosmology and σ₈ = 0.9 and another which we refer to as the ‘EdS’ simulation which has Ω_m = 1. The EdS simulation is not a viable cosmological model for our Universe as it has already been ruled out by many observations but we use it here as a tool to examine how robust our void models are to large changes in the power spectrum or cosmology.

The simulation details are summarized in Table 2. Most of the simulations use N = 256³ particles to represent the dark matter while for the larger simulation box of 500 h⁻¹ Mpc we use 400³ particles. The error on the abundance of voids measured in the 500 h⁻¹ Mpc box is estimated from eight lower resolution simulations which have 256³ particles in a computational box of 500 h⁻¹ Mpc on a side. These lower resolution simulations have a mean abundance which agrees with the 400³ particle simulation over the range of scales which we consider and are computationally less expensive to run and analyse with the void finder. The initial conditions of the particle load were set up with a glass configuration of particles (Baugh, Gaztanaga & Efstathiou 1995) and the Zel'dovich approximation to displace the particles from their initial positions. We chose a starting redshift of z = 100 in order to limit the discreteness effects of the initial displacement scheme (Smith et al. 2003). The linear theory power spectrum used to generate the initial conditions was created using the camb package of Lewis & Bridle (2002). Snapshot outputs of the dark matter distribution as well as the group catalogues were made at redshifts 1, 0.5 and 0. In the following section, we also test voids that are identified with dark matter haloes. The simulation code gadget-2 has an inbuilt friends-of-friends (FOF) halo finder which was applied to produce halo catalogues of dark matter particles with 10 or more particles. A linking length of 0.2 times the mean interparticle separation was used in the halo finder.

Table 2.

Details of the simulations used in this work.

Model	Ωm	h	σ8	Box size	No. of particles	No. of simulations	z output
ΛCDM	0.26	0.715	0.8	64, 128, 256	2563	8	0, 0.5, 1
				500	4003	1	0
				500	2563	8	0
ΛCDM	0.26	0.715	0.9	64, 128, 256	2563	8	0, 0.5, 1
EdS	1	0.7	0.8	64, 128, 256	2563	8	0

4 RESULTS

In the following sections, we compare simulation results for the abundance of voids with the predictions of both the volume-conserving and the SVdW model. In Section 4.1, we present the measured abundances of voids from the ΛCDM simulations in four different simulation box sizes. In Sections 4.2–4.4, we test the robustness of the models to variation in the critical void underdensity, redshift and cosmology, respectively. In Section 4.5, we present the abundance of voids identified in the dark matter halo catalogue.

4.1 Baseline model comparison

We implement the void finder, which is described in Section 3.1, to measure the abundance of spherical voids which have ρ_v = 0.2ρ_m at z = 0 in all four simulation box sizes of the ΛCDM cosmology, see Table 2. Fig. 6 shows the average number density as a function of radius, of voids measured from eight different realizations of the ΛCDM cosmology in simulation box sizes 64 h⁻¹ Mpc (green), 128 h⁻¹ Mpc (purple), 256 h⁻¹ Mpc (red) and 500 h⁻¹ Mpc (cyan) on a side. The error bars represent the scatter amongst these simulations.

Figure 6.

Void abundance in simulations versus predictions with ρ_v = 0.2ρ_m in the dark matter distribution of the σ₈ = 0.8 ΛCDM cosmology in simulation box sizes 64 h⁻¹ Mpc (green), 128 h⁻¹ Mpc (purple), 256 h⁻¹ Mpc (red) and 500 h⁻¹ Mpc (cyan) on a side. The error bars represent the scatter on the mean from eight different realizations of this cosmology in each box size. The range in predictions cover the parameter interval δ_c = [1.06,1.686] with δ_v = −2.7 and are consistent with simulations for V dn (grey shaded) but not SVdW models (orange hatched).

The orange hatched region in this figure represents the SVdW model within the parameter interval δ_c = [1.06,1.686] and δ_v = −2.7 and assuming that the voids have expanded by a factor of 1.7 today. The grey shaded region shows the V dn model for the same parameters. As discussed in Section 2.2, the range in δ_c accounts for the void-in-cloud process by which a void in a larger overdense regions will be crushed out of existence. As we can see from Fig. 6, this only affects the smallest voids of r < 1 h⁻¹ Mpc and for larger voids the abundance is insensitive to δ_c. The decrease in the void abundance at r(h⁻¹Mpc) ∼ 2.5, 1.5 and 1 for the 256, 128 and 64 h⁻¹ Mpc boxes shows the resolution limit for each of these simulations where small voids are not fully resolved and so the abundance is decreased.

As a result of assuming an isolated spherical expansion model, the SVdW model overpredicts the abundance by a factor of 5 whereas the V dn model agrees with simulations to ∼16 per cent across the range 1 < r(h⁻¹ Mpc) < 15 where the results measured from simulations in different box sizes have converged. This shows that the excursion set model is in good agreement with simulations once we account for the fact that voids merge as they expand and do not conserve the linear theory number density. The V dn model conserves the volume rather than the number of voids and hence implies that the number density decreases in going from the linear to the non-linear regime by the same amount that the volume of the voids grow. It is somewhat surprising that using the factor of 1.7 in this model, which applies to the expansion of isolated objects, fits the results from the simulations where voids have merged as they expand. It is important to test that this is not just a coincidence but rather is robust to other choices of parameters in the simulations.

4.2 Underdensity variation

In both the SVdW and V dn models, we adopt the shell-crossing criterion ρ_v = 0.2ρ_m for defining the void and match predictions to ρ_v as defined by the simple void finder of Section 3.1. If the agreement between the V dn model and simulations was robust, we would expect that it would be preserved for at least small variations in this criterion.

We modify our void finder such that the largest non-overlapping spherical regions which have densities ρ_v = 0.3ρ_m and ρ_v = 0.4ρ_m are recovered from the simulations. The results are shown in the left- and right-hand panels of Fig. 7, respectively. The errors plotted in this figure represent the scatter on the mean from eight simulations.

Figure 7.

Void abundance for different defining underdensities ρ_v = 0.3ρ_m (left-hand panel), ρ_v = 0.4ρ_m (right-hand panel). The grey shaded region represents the excursion set predictions with varying amplitude and using a linear underdensity value δ_c = 1.686 and δ_v, given in the legend in each panel. The amplitude rescaling versus the SVdW predictions ranges from ρ_v/ρ_m (V dn; top black dashed curve) to 1/5 (bottom black dotted curve) both of which preserve agreement for ρ_v = 0.2ρ_m.

As discussed in Section 2.2 (see also Fig. A1 and equation B1), changing the underdensity criteria in the spherical evolution model alters the shape of the abundance function through the linear threshold δ_v. Specifically, for ρ_v = 0.3ρ_m, δ_v = −1.8; while for ρ_v = 0.4ρ_m, δ_v = −1.24.

For dark matter voids with ρ_v = 0.2ρ_m, the predicted abundance for the V dn model are approximately a factor of 5 smaller than those of the SVdW model. In modelling the number density of underdense regions with ρ_v = 0.3ρ_m and ρ_v = 0.4ρ_m, which cannot be directly related to shell crossing in the spherical expansion model, we adopt a phenomenological approach. Given that for ρ_v = 0.2ρ_m, the V dn model has the same shape as the SVdW model but a factor of ρ_v/ρ_m = 1/5 lower amplitude, we can preserve the good match there either by following the V dn prescription literally and rescaling SVdW by ρ_v/ρ_m or by simply keeping this factor fixed at 1/5.

This range is plotted in Fig. 7 as grey shaded regions bounded by the two limiting cases, black dashed and dotted lines. Simulation results clearly favour the simple phenomenological prescription of rescaling the amplitude by 1/5. The ρ_v/ρ_m scaling prescription of the V dn model would overpredict the amplitude by ∼1.5 for voids with ρ_v = 0.3ρ_m and 2 for ρ_v = 0.4ρ_m. These results again highlight the point that excursion set models predict the overall shape of abundance function accurately and only the amplitude needs to be altered to fit the simulations results, here without the benefit of volume conservation as motivation. Note that the preferred rescaling of 1/5 is more than sufficient to bring the predictions to a physical void filling fraction for ρ_v ≥ 0.2ρ_m.

4.3 Redshift variation

Next, we check the robustness of results to the redshift at which the void abundance is measured. Fig. 8 shows the number density of voids as a function of radius at z = 0.5 (blue) and z = 1 (red) measured from the ΛCDM, σ₈ = 0.8 simulation. The measured abundances from the three simulation box sizes 64, 128 and 256 h⁻¹ Mpc are the volume-weighted averages and errors over eight realizations. The volume-conserving (V dn) model is shown as a black hatched (grey shaded) region using δ_v = −2.7 at z = 0.5 (z = 1) and the parameter range δ_c = [1.06, 1.686]. Note that we have used only one colour for the results from the three simulation boxes at each redshift for clarity in this figure.

Figure 8.

Redshift dependence of the void abundance with ρ_v = 0.2ρ_m at z = 0.5 (blue) and z = 1 (red) measured from the ΛCDM, σ₈ = 0.8 simulation. The black hatched (grey shaded) region represents the V dn model using linear underdensity values of δ_v = −2.7 at z = 0.5 (z = 1) for the range δ_c = [1.06, 1.686]. Note that the measured average abundances and errors from the 256, 128 and 64 h⁻¹ Mpc simulation boxes are the volume-weighted values. Note also that in this figure we have plotted the results from the three simulation boxes using the same colour for clarity.

We again find that the V dn model works very well in reproducing the abundance of voids in the dark matter in a ΛCDM universe at both redshifts, while the SVdW model, which is not plotted here for clarity, again overpredicts the abundances by approximately a factor of 5. Fig. 8 shows that smaller (larger) voids are more (less) abundant at z = 1 compared to z = 0.5 which is also found in the model predictions at both redshifts.

4.4 Cosmological parameter variation

In order to check if the V dn model for the abundance of voids works when we change the cosmological model, we have run two simulations of alternative cosmologies to the standard ΛCDM with σ₈ = 0.8 which we discussed in the previous section. In the first alternative cosmology, we have chosen to modify only the value of σ₈ to 0.9, see Table 2; our second alternative cosmology is an EdS universe where the matter density parameter Ω_m = 1. The linear perturbation theory power spectra for these simulations were generated using camb (Lewis & Bridle 2002) and normalized to σ₈ = 0.9 (σ₈ = 0.8) for the ΛCDM (EdS) simulations in order to generate the initial conditions for the simulations and the variance σ(R) which is used in the excursion set model for the abundance.

The measured z = 0 number density of voids with ρ_v = 0.2ρ_m in the σ₈ = 0.9 and EdS simulations are shown in Fig. 9. The volume-conserving model is shown in both panels as a grey shaded region as in previous plots. We have used the same value, δ_v = −2.7, for the linear perturbation theory underdensity. Note that this parameter is different in different cosmologies; however, we find that such a small change in δ_v going from an EdS to a ΛCDM universe has a small impact on the predicted abundance of voids in the excursion set theory and the main differences arise from the change in the variance, σ(R) (see Appendix A).

Figure 9.

Void abundance for alternate cosmological parameters at z = 0. Upper: ΛCDM with the initial conditions normalized to give σ₈ = 0.9. Lower: EdS model with σ₈ = 0.8 with Ω_m = 1. The grey shaded region shows the V dn model within the parameter interval δ_c = [1.06,1.686] and using δ_v = −2.7.

From Fig. 9, it is clear that the volume-conserving model works well in both of these cosmologies and fits the abundance of voids to within 25 per cent over the range 1 < r(h⁻¹ Mpc) < 15. It is interesting to note the overall decrease in the abundance of voids in the dark matter distribution for voids with small radii r < 2h⁻¹ Mpc in these two cosmologies which is most obvious in the measured number density from the EdS simulation and a larger abundance for the σ₈ = 0.9 cosmology for large r. It is also clear from Fig. 9 (lower) that the excursion set model predicts more squashing of smaller voids due to the void-in-cloud effect but this is occurring right on the resolution limit of our simulations at 2 < r(h⁻¹ Mpc). Finally, note that even if we modified the SVdW model in the ad hoc manner of Fig. 4 to match the simulation results of ΛCDM with σ₈ = 0.8, the predictions would be far off simulation results for the EdS cosmology.

4.5 Halo defined voids

Voids in the galaxy population are defined not through the dark matter density field but by the number density field n_h of the dark matter haloes they populate. In this section, we use density minima in the halo number density. Our goal is to test how faithfully the abundance of voids in the dark matter matches that in the halo populations within the context of the simple void finder of Section 3.1. It is important to note that a comparison between voids in the dark matter and halo distributions should account for the galaxy/dark matter halo biasing relation. Benson et al. (2003) showed that the properties of galaxy and dark matter voids differ significantly as a results of galaxy bias e.g. if galaxies are sparse tracers of the underlying dark matter, this gives rise to larger voids in the galaxy distribution. Furlanetto & Piran (2006) also studied the abundance of voids in the galaxy distribution within the excursion set formalism and showed that after accounting for bias, galaxy voids should be larger than dark matter voids, while voids selected using luminous galaxies should be larger than those using faint galaxies.

In this section, we use the FOF halo catalogues from the 128, 256, 500 h⁻¹ Mpc simulation boxes and the publicly available halo catalogues from the MultiDark and Bolshoi simulations (Riebe et al. 2011) which have computational box sizes of L = 1000 and 250 h⁻¹ Mpc on a side, respectively. These haloes have been identified using the bound density maxima algorithm (Klypin & Holtzman 1997). We only use haloes which have V_max > 200 km s⁻¹ and M > 1 − 2 × 10¹² h^{− 1} M_⊙ from the Bolshoi and MultiDark simulations to ensure that the statistics are robust. We use the void finder described in Section 3.1 to identify voids in the distribution of haloes which have n_v = 0.2n_h where n_v is the average number density in the void whereas n_h is the average in the whole simulation. Our final sample consists of 5768 voids using 1.7 × 10⁶ haloes from the MultiDark simulation and 4826 voids using 2.2 × 10⁶ haloes from the Bolshoi simulation. Both of these simulations are of a higher resolution than the ones we carried out in 128, 256, 500 h⁻¹ Mpc simulation boxes – it is useful to compare the abundance of voids in the halo population from these simulations to ours as an indication of the scales at which our results have converged.

The measured abundance of voids in the haloes population from our simulations and the Bolshoi and MultiDark catalogues at z = 0 are shown in Fig. 10. The errors shown on the results from the 128 (purple), 256 (red), 500 (cyan) h⁻¹ Mpc simulation boxes are measured from the scatter amongst eight different realisations in each box size. The errors on the MultiDark simulation were obtained by Jackknife sampling from each simulation by dividing the simulation volume into N_sub = 8 equal subvolumes and then systematically omitting one subvolume at a time in order to calculate the void abundance on the remaining N_sub − 1 volume. We find that voids identified in this manner through the halo distribution do not follow the V dn model assuming n_v = 0.2n_h which corresponds to dark matter voids of ρ_v = 0.2ρ_m. They also do not follow the SVdW model which would have the same shape but five times the amplitude.

Figure 10.

The number density of voids with n_v = 0.2n_h in the halo distribution from the 128 (purple), 256 (red), 500 (cyan) h⁻¹ Mpc simulation boxes. The results from the MultiDark (Bolshoi) simulation are shown in dark green (blue). The error bars represent the error on the mean from eight simulations. The errors on the MultiDark simulation represent the Jackknife error on the mean. The grey shaded region bounded by the black dashed and dotted line represents the volume-conserving model with δ_v = −1.24 and varying amplitude as in Fig. 7. The grey solid line represents the V dn model with δ_v = −2.7.

Overall Fig. 10 shows that our void finder finds large halo defined voids that do not correspond to dark matter defined voids of the same underdensity for r ≳ 10 h^{− 1} Mpc. Although it is difficult to compare these results with previous work due to the large differences in the void finders used, qualitatively this agrees with the findings of Benson et al. (2003) who measured the VPF from simulations and found that the VPF for voids with r > 5 h⁻¹ Mpc was much higher for the galaxy catalogues compared to the dark matter. These results illustrate the fact that there is not always a 1:1 correspondence between voids in the dark matter and the dark matter halo distributions and this is especially pronounced when we define a void as having a fixed underdensity which is the same for dark matter and haloes.

To illustrate the mismatch between the voids which we find in the dark matter and halo distributions using the same underdensity criterion, Fig. 11 shows a 10 × 50 × 50 h⁻¹ Mpc (left-hand panel) and a 60 × 14 × 60 h⁻¹ Mpc (right-hand panel) slice through the dark matter density field which has been evaluated on a grid of 256³ points from the 500 h⁻¹ Mpc simulation box. The coloured contours represent the log of the density field in each cell and the haloes around each void are represented by black dots. The radius of each void is r ∼ 21 h⁻¹ Mpc (left-hand panel) and r ∼ 26 h⁻¹ Mpc (right-hand panel) and is shown as a grey dashed line in Fig. 11. The red circles in this plot show the voids identified in the dark matter whose centres are within 10 h⁻¹ Mpc of the centre of the void in the halo population. Not only is it possible to find more than one dark matter void which overlaps with the halo void but the radii of the dark matter voids at which ρ_v = 0.2ρ_m are a lot smaller than the halo voids which satisfy the analogous criterion.

Figure 11.

Left: a 10 × 50 × 50 h⁻¹ Mpc slice through the 500 h⁻¹ Mpc simulation box centred on a large, r ∼ 21h⁻¹ Mpc, void in the halo population (black dots). The diameter of the void is shown as a dashed grey line and the coloured contours represent the log of the density field which has been evaluated on a grid of 256³ points. The red circles represent all the voids in dark matter which have ρ_v = 0.2ρ_m and whose centres are within 10 h⁻¹ Mpc of the centre of the void in the halo distribution. Right: a 60 × 14 × 60 h⁻¹ Mpc slice through the 500 h⁻¹ Mpc simulation box centred on a large, r ∼ 26h⁻¹ Mpc, void in the halo population (black dots). The red circles represent all the voids in dark matter which have ρ_v = 0.2ρ_m and whose centres are within 10 h⁻¹ Mpc of the centre of the void in the halo distribution. Note that these voids only appear to be overlapping due to the projection effect.

There are at least two possible ways to reconcile the V dn predictions for the abundance of dark matter voids with that of the halo voids. First, a scale-dependent modification to the barrier in the V dn model could be used to alter the underdensity threshold used to find voids in the dark matter. Secondly, if we keep a fixed underdensity threshold to define dark matter voids, it may be possible to find a scaling of this threshold to define voids in the halo distribution. These ideas are beyond the scope of this work but see Furlanetto & Piran (2006) for related ideas.

As a simpler illustration of these ideas, in Fig. 10 we also plot the V dn model assuming that halo defined voids of n_v = 0.2n_h correspond to dark matter defined voids of ρ_v = 0.4ρ_m. These predictions are plotted as a grey shaded region allowing the amplitude to vary from the predictions of the V dn model which rescales the SVdW amplitude by ρ_v/ρ_m = 0.4 (black dashed lines) and the rescaling of 1/5 that fits our dark matter voids well as in Fig. 7 (black dotted line). Compared to the predictions of the V dn model for dark matter voids of n_v = 0.2n_h (solid grey line), these black dashed and dotted curves match the abundance of voids in the halo populations better though no single rescaling matches perfectly across the full range.

5 SUMMARY AND CONCLUSIONS

The next generation of galaxy redshift surveys such as BigBOSS (Schlegel et al. 2009), Euclid (Laureijs et al. 2011) and WFIRST (Albrecht et al. 2006; Green et al. 2012) will allow us to study the large-scale structure of our Universe in ever greater detail. Cosmic voids represent one of the main components which strongly influence the growth of clusters, walls and filaments in the mass distribution. Studying the statistics and dynamics of these underdense regions is a promising way to test the cosmological model and hierarchical structure formation.

Several challenges which may affect the usefulness of voids as a probe of cosmology are addressed in this paper, such as, given a survey or numerical simulation of a given size, how robust are the statistics on the number density of voids of a given size, how accurately can we predict the number density of these voids and how faithfully do voids in the halo population represent voids in the dark matter. Using N-body simulations of a ΛCDM cosmology, we test the excursion set model for the abundance of voids including the model provided by Sheth & van de Weygaert (2004), which takes into account the void-in-void and void-in-cloud scenarios. Our void finder makes use of the zobov (Neyrinck 2008) algorithm which uses Voronoi tessellation to locate density minima. We define a void as a spherical region around these minima with ρ_v = 0.2ρ_m and make use of several different computational box sizes so we can determine the volume and resolution that is needed in order to recover robust statistics for voids of a given size. We have tested the robustness of our void finder to the following changes and have found convergent results: using different simulation box sizes and particle numbers, using a volume-weighted centre or the core particle to define the centre of a void, using all particles around the density minima or only particles in a zone given by zobov and using only statistically significant voids or voids with a core particle density <0.2ρ_m.

We find that the measured abundance of voids at z = 0 from a ΛCDM simulation does not match the predictions of the Sheth & van de Weygaert (2004) model, which greatly overpredicts the results using a linear underdensity of δ_v = −2.7. We find that the excursion set theory, which is the basis of the SVdW model, accurately predicts the shape of the abundance of voids measured from simulations. However, the predicted amplitude is incorrect due to the assumption of isolated spherical expansion which does not account for the merging of voids as they expand. Instead, we find a volume-conserving model, which is also based on the excursion set method with δ_v = −2.7, matches the measured abundances to within 16 per cent for void radii 1 < r(h⁻¹ Mpc) < 15. This model works remarkably well and suggests that the number density of voids decreases in going from the linear to the non-linear regime by the same amount that the voids expand. This agreement is robust to varying the redshift in the ΛCDM model as well as the underlying cosmology. Using simulations of different cosmological models, a ΛCDM cosmology with σ₈ = 0.9 and a EdS cosmology with Ω_m = 1, we find that the volume-conserving model works well and reproduces the measured number density from each simulation to within 25 per cent over the range 1 < r(h⁻¹ Mpc) < 15. We have also tested model predictions for density thresholds of ρ_v = 0.3ρ_m and ρ_v = 0.4ρ_m and find that the volume fraction physicality rescaling factor remains fixed at ∼1/5 rather than scaling as ρ_v/ρ_m.

Using the number density threshold criterion of n_v = 0.2n_h in our void finder, we have examined the voids in the halo population from the 128, 256 and the 500 h⁻¹ Mpc on a side computational boxes. We also use the Bolshoi and the MultiDark (Riebe et al. 2011) simulations to measure void abundances. These two simulations are of a higher resolution than our simulations and we have confirmed that our measured void abundances in the halo population agree with both the MultiDark and Bolshoi measurements. We find that the number density of voids in the dark matter halo distribution is very different to that in the dark matter, with fewer small, r < 10 h⁻¹ Mpc, voids and many more large, r > 10 h⁻¹ Mpc, voids. These results indicate that a given void in the halo distribution of fixed underdensity of n_v = 0.2n_h cannot be unambiguously related to a void in the dark matter of equal underdensity and in the case that there is a 1:1 correspondence, the radii at which a dark matter or halo void have a given underdensity can be very different.

Cosmic voids are a promising and interesting tool which can be used to test many aspects of the ΛCDM cosmological model and having an accurate model for the number density of voids in the Universe represents a first step. In this work, we have presented a model based on the excursion set theory which conserves the volume fraction of voids and works well at predicting the abundance of voids identified in dark matter from N-body simulations over a range of scales. Establishing the relationship between voids in the dark matter and in the halo or galaxy distribution requires further study.

EJ acknowledges the support of a grant from the Simons Foundation, award number 184549. This work was supported in part by the Kavli Institute for Cosmological Physics at the University of Chicago through grants NSF PHY-0114422 and NSF PHY-0551142 and an endowment from the Kavli Foundation and its founder Fred Kavli. YL and WH were additionally supported by US Deptartment of Energy contract DE-FG02-90ER-40560 and the David and Lucile Packard Foundation. We are grateful for the support of the University of Chicago Research Computing Center for assistance with the calculations carried out in this work.

APPENDIX A: Spherical evolution model

In this appendix, we review the spherical evolution model which describes the non-linear evolution of an un-compensated spherically symmetric tophat underdense (overdense) perturbation. To illustrate the physics, we will use the spherical tophat model in EdS cosmology as an example, which is analytically solvable until shell crossing.

Consider an initial spherical tophat density perturbation (|δ₀| ≪ 1) of physical radius R₀ at a_i = a(t_i). We can think of the perturbation as composed of concentric mass shells, labelled by their initial physical radii R_i at a_i. Let Δ(R_i, a) denote the average overdensity of the region enclosed by the mass shell R_i. Then, its initial value is

\begin{eqnarray} \Delta _{\rm i}(R_{\rm i}) \equiv \Delta (R_{\rm i}, a_{\rm i}) = \left\lbrace \begin{array}{ll}\delta _0 & \quad R_{\rm i} \le R_0, \\ \delta _0 (R_0/R_{\rm i})^3 & \quad R_{\rm i} > R_0. \end{array} \right. \end{eqnarray}

(A1)

For brevity, we omit the R_i argument of Δ_i below. According to Birkhoff's theorem, the evolution of the mass shell R_i only depends on the total mass inside R_i, but not the mass distribution or the mass outside. Thus, the shell R_i evolves in the same way as an FLRW universe

\begin{eqnarray} \left[\frac{\dot{R}(t;R_{\rm i})}{R(t;R_{\rm i})}\right]^2 = H_{\rm i}^2 \left[(1 + \Delta _{\rm i}) \left( \frac{R_{\rm i}}{R} \right)^{3} - \frac{5}{3} \Delta _{\rm i} \left( \frac{R_{\rm i}}{R} \right)^{2} \right], \end{eqnarray}

(A2)

where R(t) is the physical radius, and initial conditions are set to the growing mode in linear theory. By introducing the dimensionless conformal time

\begin{eqnarray} {\rm d}\eta = \frac{R_{\rm i}}{R} \sqrt{ \left| \frac{5}{3} \Delta _{\rm i} \right| } \, H_{\rm i} {\rm d}t \,, \end{eqnarray}

(A3)

we can solve equation (A2) in a parametric form (to leading order in Δ_i)

\begin{eqnarray} \frac{R}{R_{\rm i}} & \simeq & \frac{1}{2} \left|\frac{5}{3} \Delta _{\rm i} \right|^{-1} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(\cosh \eta - 1) &\quad \delta _0<0, \\ (1 - \cos \eta ) &\quad \delta _0>0; \\ \end{array}\right. \end{eqnarray}

(A4)

\begin{eqnarray} H_{\rm i} t & \simeq & \frac{1}{2} \left| \frac{5}{3} \Delta _{\rm i} \right|^{-3/2} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(\sinh \eta - \eta ) &\quad \delta _0<0, \\ (\eta - \sin \eta ) &\quad \delta _0>0.\\ \end{array}\right. \end{eqnarray}

(A5)

Given this evolution of mass shells, we are particularly interested in shell crossing for the expansion case (δ₀ < 0), which is usually seen as a characteristic event that signifies the formation of the void at a non-linear level. Note that these solutions represent a family of trajectories labelled by R_i and parametrized by |$\eta _{R_{{\rm i}}}$|⁠. We can find out when and where shell crossing first occurs by differentiating the parametrized solutions with respect to R_i and η, and requiring that dR and dt vanish, for R_i > R₀

\begin{eqnarray} {\left[\begin{array}{cc}A_{11} &\quad A_{12} \\ A_{21} &\quad A_{22} \end{array}\right]} {\left[\begin{array}{cc}{\rm d}R_{\rm i}/R_{\rm i} \\ {\rm d}\eta \end{array}\right]} = 0, \end{eqnarray}

(A6)

where

\begin{eqnarray} A_{11} &=& 2 \left| \frac{5}{3} \Delta _{\rm i} \right|^{-1} (\cosh \eta - 1),\nonumber \\ A_{12} &=& \frac{1}{2} \left| \frac{5}{3} \Delta _{\rm i} \right|^{-1} \sinh \eta , \nonumber \\ A_{21} &=& \frac{9}{4} \left| \frac{5}{3} \Delta _{\rm i} \right|^{-3/2} (\sinh \eta - \eta ), \nonumber \\ A_{22} &=& \frac{1}{2} \left| \frac{5}{3} \Delta _{\rm i} \right|^{-3/2} (\cosh \eta - 1). \end{eqnarray}

(A7)

For this homogeneous system of linear equations to have non-zero solutions, we must have |$\det A=0$|⁠. Thus, we derive the shell-crossing condition

\begin{eqnarray} \frac{\sinh \eta \, (\sinh \eta -\eta )}{(\cosh \eta -1)^2} = \frac{8}{9}. \end{eqnarray}

(A8)

Shell crossing first happens at η_sc = 3.488 among the boundary shells, i.e. R_i = R₀ in the above criterion. At shell crossing, the void interior has a relative density

\begin{eqnarray} 1+\Delta _{\rm sc} \simeq \frac{9}{2} \frac{(\sinh \eta _{\rm sc}-\eta _{\rm sc})^2}{(\cosh \eta _{\rm sc}-1)^3} = 0.2047 \,, \end{eqnarray}

(A9)

which implies that the void has expanded by a factor of (1 + Δ_sc)^{− 1/3} = 1.697 in comoving radius. Note that these numbers do not depend on the size of the void.

To calculate the linear theory prediction of the void underdensity at shell crossing δ_v, we expand R(t) to the first order with the help of the parametric solution (A4) and (A5), for R_i ≤ R₀

\begin{eqnarray} \frac{R}{R_{\rm i}} = \frac{a}{a_{\rm i}} \left[1 - \frac{\delta _0}{3} \left(\frac{3}{2} H_{\rm i} t\right)^{-2/3} + \cdots \right], \end{eqnarray}

(A10)

the first order of which gives the linear underdensity

\begin{eqnarray} \delta = \delta _0 \left(\frac{3}{2} H_{\rm i} t\right)^{-2/3} \simeq - \frac{3}{20} \big [6(\sinh \eta -\eta )\big ]^{2/3}. \end{eqnarray}

(A11)

Thus, the linear underdensity at shell crossing is

\begin{eqnarray} \delta _{\rm v} = \delta (\eta _{\rm sc}) = -2.717 \,. \end{eqnarray}

(A12)

Note that this number is different in different cosmologies, e.g it is δ_v = −2.731 for ΛCDM. However, we find that such a small change in δ_v going from an EdS to a ΛCDM universe has a small impact on the predicted abundance of voids in the excursion set theory. Also note the value for δ_v for the EdS universe quoted here is a correction to the result presented in Sheth & van de Weygaert (2004) (δ_v = −2.81).

Similarly, for the spherical collapse model

\begin{eqnarray} \delta = \delta _0 \left(\frac{3}{2} H_{\rm i} t\right)^{2/3} \simeq \frac{3}{20} \big [6(\eta -\sin \eta )\big ]^{2/3}. \end{eqnarray}

(A13)

The well-known turn-around and virialization of haloes occur at η_ta = π and η_vir = 2π, leading to δ_c = 1.062 and δ_c = 1.686, respectively. For the ΛCDM model, these become δ_c = 1.303 and δ_c = 1.674, for turn-around and collapse at z = 0. The EdS range encompasses that of the ΛCDM parameters and so in the main text we have adopted the EdS parameters to show the full range of possibilities.

Finally, note that these linear density thresholds δ_v and δ_c, which are to be used in the excursion set formalism, are independent of the size of the structures.

APPENDIX B: SVdW model modifications

In this appendix, we explore modifications of the SVdW model prescription and approximations introduced in Sheth & van de Weygaert (2004). The spherical evolution model relates the linear underdensity δ to non-linear underdensity Δ, or alternatively to r/r_L = (1 + Δ)^−1/3, where r is the void radius r = R(t; R₀)/a(t) and r_L the linear radius r_L = R₀/a_i. If we relax the criterion for defining a void to correspond to underdense regions that have undergone shell crossing, there is additional freedom in defining the void abundance as a function of radius so long as δ_v and r/r_L are chosen self-consistently.

We show this relation before shell crossing for the EdS model in Fig. A1. The relation is well fitted to (Bernardeau 1994)

\begin{eqnarray} \delta _{\rm v}\approx c \left[1 - (r/r_{\rm L})^{3/c}\right] = c \left[1 - (\rho _{\rm v}/\rho _{\rm m})^{-1/c}\right], \end{eqnarray}

(B1)

where c = 1.594, with errors below 0.2 per cent. Also shown is the maximum r/r_L, constrained by requiring the total volume fraction in voids from equation (14) be less than 1, in the SVdW model. Clearly, this constraint depends on δ_c, which in the plot is chosen to be the least stringent value 1.06 in the expected range. Note that no choice of δ_v and r/r_L is physical for they all violate the total volume condition. In the main text, we also considered ad hoc modifications of the model where δ_v and r/r_L are considered unrelated.

Figure A1.

Linear underdensity as a function of void expansion factor (exact: dot–dashed line, fit of equation (B1): blue solid line). Also shown is the physicality constraint from requiring that the total volume fraction in voids from equation (14) remain less than unity in SVdW model for the least stringent value of δ_c = 1.06 (shaded region). Only as |δ_v| → 0, where the regions are at the mean density and do not expand, does the model regain physicality. Bottom panel shows fractional error in the fit.

Sheth & van de Weygaert (2004) also utilized an approximation to the exact prediction for the abundance function of equation (6) which introduces notable errors for scales where the void-in-cloud process dominates and consequently the total volume fraction. Their approximation,

\begin{eqnarray} f_{\ln \sigma }(\sigma ) &\approx & \sqrt{\frac{2\nu }{\pi }}\exp \left(-\frac{\nu }{2}\right) \exp \left( -\frac{|\delta _{\rm v}|}{\delta _{\rm c}}\frac{\mathcal {D}^2}{4\nu } - 2\frac{\mathcal {D}^4}{\nu ^2}\right), \end{eqnarray}

(B2)

where |$\nu = \delta _{\rm v}^2/\sigma ^2(M)$|⁠, had a stated realm of validity of δ_c/|δ_v| > 1/4 or |$\mathcal {D}<4/5$|⁠. Unfortunately, this approximation has uncontrolled errors at ν ≪ 1, exactly where the void-in-cloud process operates as shown in Fig. B1. Our piecewise approximation is accurate at the 0.2 per cent level or better everywhere. Note that the errors and smoothness of our approximation can be improved at the transition point by a suitable interpolation between the two piecewise curves though it is not necessary for this work.

Figure B1.

SVdW model void abundance is shown with three different f_ln σ: the exact formula (6) (thick black dot–dashed line), the Sheth & van de Weygaert (2004) SVdW approximation (B2) (red dashed line) and our piecewise approximation (8) (blue solid line), along with the absolute value of the error relative to exact. We use δ_v = −2.7 and δ_c = 1.06, well within the stated domain of validity of both approximations. The SVdW approximation breaks down where the void-in-cloud process becomes important, whereas the piecewise approximation is accurate everywhere with errors peaking at 0.2 per cent where the piecewise transition occurs. We use the σ₈ = 0.8 ΛCDM cosmology as listed in Table 2.