The Santiago–Harvard–Edinburgh–Durham void comparison – I. SHEDding light on chameleon gravity tests

Abstract

1 INTRODUCTION

After billions of years of evolution since the big bang, the Universe has developed rich large-scale structures, the so-called cosmic web. This forms a complex and space-filling pattern, and is composed of lumpy knots at the intersection of long filamentary bridges, which in turn form at the intersection of tenuous sheets (Bond, Kofman & Pogosyan 1996, see Cautun, van de Weygaert & Jones 2013 for a visualization). These structures are populated by dark matter haloes, some of which are illuminated by galaxies such as our own Milky Way.

Such rich structures are a consequence of the interplay of various ingredients, most of which have their roots in fundamental physics. The tiny initial density fluctuations are believed to be seeded by quantum effects in a hypothetical period of inflationary expansion. The matter species are predicted by the standard model of particle physics, along with its extensions that provide theoretical candidates for the dark matter. The effect of gravity is felt on almost all scales, from determining the expansion rate of the Universe as a whole to governing the formation and evolution of the cosmic web itself and the galaxies and stars inside them. Finally, there is the mysterious cosmic acceleration – the observation that the expansion rate of the Universe has been accelerating instead of decelerating as predicted originally – discovered initially almost two decades ago (Riess et al. 1998; Perlmutter et al. 1999), which could be due to a small positive cosmological constant, a new exotic matter species (dark energy) or a breakdown of Einsteinian general relativity (GR) on cosmic scales (modified gravity). This interplay suggests that the accurate determination of the properties of the cosmic web could help to improve our understanding of the laws of fundamental physics that have shaped the Universe. However, it also poses great challenges to find optimal ways of extracting accurate information from a highly entangled observational reality.

One of the key properties of the cosmic web, which has been a subject of growing interest in recent years, is related to the vast space in between sheets and filaments, nearly empty of galaxies and dark matter. These are commonly referred to as ‘cosmic voids’ (van de Weygaert & Platen 2011). By definition, these are low-density regions that contain little matter and few galaxies; however, they still encode detailed information about the underlying cosmological model (e.g. Li 2011; Bos et al. 2012; Cai, Padilla & Li 2015; Massara et al. 2015; Pisani et al. 2015; Zivick et al. 2015; Achitouv 2016; Banerjee & Dalal 2016; Demchenko et al. 2016). Voids make an attractive cosmological probe since their properties are hardly affected by galaxy formation physics, which is still a major unknown, and are well reproduced by dark matter only simulations (after suitably populating haloes with galaxies) (Paillas et al. 2017). Their abundance, for example, can be qualitatively understood using semi-analytical models (see e.g. Sheth & van de Weygaert 2004; Paranjape, Lam & Sheth 2012; Jennings, Li & Hu 2013). However, as we will discuss below, there is an uncertainty in the precise definition of voids, which inevitably makes void abundance a less well-defined quantity. In recent years, other observables associated with voids have attracted attention, e.g. void redshift-space distortions, which is the anisotropy in the observed correlation of voids and galaxies in redshift space induced by the peculiar velocities of galaxies (Hamaus et al. 2015, 2016; Cai et al. 2016). An additional observable is void lensing, an antilensing effect caused by the underdense voids deflecting light away from their centres, which has been detected in both the Sloan Digital Sky Survey (sdss) and the Dark Energy Survey (des) (Clampitt & Jain 2015; Gruen et al. 2016; Sanchez et al. 2017).

The surge of interest in cosmic voids is also partly thanks to a community that works on modified gravity theories (see e.g. Joyce et al. 2015 and Koyama 2016 for the latest reviews). This growing field explores the possibility that the cosmic acceleration is due to new gravitational physics on cosmological scales and attempt to use precision cosmological data to verify the validity of GR in a new regime, beyond traditional tests of the theory (Will 2014). In environments similar to the Solar system (e.g. deep or large gradient of the gravitational potential, high matter density), small-scale tests of GR have placed stringent constraints on the viability of any new gravity theory (or extension to GR). As a result, many of the newly proposed theories have a so-called screening mechanism, by which standard GR is recovered, and existing small-scale tests of GR satisfied. As a result of this, the differences between these models and GR are usually suppressed in dense regions where most cosmological data are currently available. It is because of this that voids offer an alternative and promising venue to test them. In the literature, there have been many studies trying to understand the properties of voids in non-standard cosmological models. While simple semi-analytical models (see e.g. Clampitt, Cai & Li 2013; Lam et al. 2015, ) can offer useful insight into the essential physics, most studies so far have relied on numerical simulations to get more precise predictions for a larger range of void properties (e.g. Li 2011; Bos et al. 2012; Li, Zhao & Koyama 2012b; Barreira et al. 2015; Cai, Padilla & Li 2015; Pisani et al. 2015; Zivick et al. 2015; Barreira et al. 2017; Paillas et al. 2017; Falck et al. 2018).

The name ‘cosmic voids’ is used to describe generic low-density regions in the Universe, and it is difficult to make a precise definition without having a clear idea about how to determine the spatial boundary of such regions. For example, looking at a three-dimensional (3D) distribution of galaxies, it is immediately apparent that regions devoid of galaxies usually take on irregular shapes and there is no straightforward way to define their spatial domains. In the literature, two frequently used methods to identify voids are based on ‘spherical underdensity’ (e.g. Padilla, Ceccarelli & Lambas 2005) and ‘space tessellation’ (e.g. Platen, van de Weygaert & Jones 2007; Neyrinck 2008). In the former approach, one centres voids at positions of lowest density and grows spheres until some density threshold, ρ_v is attained (e.g. the average density inside the sphere is less than ρ_v). In the second approach, one covers the space with tessellation cells (such as Voronoi and Delaunay) in such a way that larger cells correspond to low-density regions. The cells are then joined to form voids according to certain prescriptions. This broad distinction, however, does not cover all possible algorithms of void finding, and actually allows for many variations in the exact implementation. For example, the threshold ρ_v in spherical underdensity void finders can take different values, the method for joining tessellation cells to form voids can vary the tracers (e.g. which galaxy population and what tracer number density) used to identify voids differ in different works, and so on. Such ambiguities in void finding could be seen as a disadvantage: if different groups cannot agree on how voids are defined, and the results from them display discrepancies among each other, how can any of these results be useful?

However, it has been suggested (see e.g. Cai et al. 2015) that a positive point of view could be taken: the cosmic web is a highly complicated structure and one should not expect that a single statistic can characterize it completely. Rather, each void-finding algorithm represents a different way to quantify the emptiness – there is no a priori way to know which algorithm is the ‘best’, since this will depend on both the observables being used and the questions being asked. The latter consideration adds a new dimension to void studies and is a main motivation of this work: given a certain cosmological model such as a specific modified gravity (which usually has distinct physical properties that imprint different observational features compared to other models), what is the void finder that best captures these differences and hence presents the best opportunity to tell it apart from other models? To address this question, in this paper we will compare the predictions of some key void properties using a number of void finders.

This is the first of a series of papers, with which we endeavour to make a detailed exploration of the performances of various void finders in predicting various void properties for various theoretical models. The emphasis here is that the different void finders should be compared (1) in specific and clear contexts, i.e. the comparison is done for certain theoretical models so that the conclusion would apply to those models and not necessarily to all models, (2) at equal footing, i.e. wherever possible, the properties of voids identified by different void finders are calculated using the same pipelines to exclude spurious differences that can bias the comparison, and (3) with a reasonable degree of reality, i.e. we focus on voids identified from galaxy fields in typical current and near-future galaxy surveys and assume that the different models being tested all predict the same galaxy clustering – the observed one. In practice, this is achieved by tuning the parameters governing how galaxies populate dark matter haloes; we do this because there is only one observed Universe, and if the galaxy clustering predicted by a model does not agree with observations, then that model should have already been ruled out. Some of the void-finding algorithms being compared here are new, and we hope that this further adds to the usefulness of this work to the community. We want to stress that it is not our intention to distinguish good versus poor void finders, but rather to find the right ones for specific questions at hand.

The only other comparison study of cosmic voids, to the best of our knowledge, is the Aspen Void Finder Comparison Project (AVFCP) by Colberg et al. (2008), which compared a total of 13 void finders; it differs from this work in several ways. First, AVFCP also compared voids identified from the dark matter field, while we only use galaxies as tracers since these are directly observable in galaxy surveys. Secondly, AVFCP identify voids from a selected region of size 60 h⁻¹ Mpc from the Millennium simulation, while we use a whole simulation box of size 1024 h⁻¹ Mpc to get a large number of voids spanning a range of sizes and study different void observables. Finally, as we have mentioned above, this is a void code comparison project in the context of testing gravity, and thus we run our void finders on several sets of simulations that differ (only) by their underlying gravity theories; in particular, we have tuned our mock galaxy catalogues so that they have the same correlation functions to avoid the contamination of having different tracers for the different models. This work will focus on the comparison of void finders in testing a particularly popular class of modified gravity models: the so-called chameleon model.

This paper is organized as follows. In Section 2, we briefly review the f(R) model, which is a representative example of chameleon models, and the simulations and mock galaxy catalogues used in the comparison. In Section 3, we introduce the six void-finding algorithms to be studied and present a visual comparison of the various void catalogues. The main results of this paper are presented in Section 4, where we compare the void abundance, void galaxy number density profile, matter density profile, and lensing tangential shear profiles; we also estimate the signal-to-noise ratio (S/N) for each void finder to determine how well they are able to discriminate this particular class of models from GR. Finally, we discuss and conclude in Section 5.

Throughout this paper, we use the unit convention c = 1, with c the speed of light, except where c is explicitly included. We use an overbar to denote the background value, and a subscript ₀ to denote the present-day value of a quantity, unless otherwise stated.

2 THEORY

We start by presenting a brief description of the f(R) gravity model studied in this paper (Section 2.1), the simulations used and the catalogues of tracers of the large-scale structure that are needed for finding voids (Section 2.2). Many of the details – such as the model and simulations – can be found elsewhere; therefore, this section is mainly for completeness and will be kept concise.

2.1 Theory, simulations, and halo/galaxy catalogues

2.1.1 f(R) gravity

In some sense, the scalar field f_R plays the role of the potential of a fifth force that is the difference between the modified gravitational force as predicted by equation (5) and the standard Newtonian force produced by the same matter density field ρ_m. There exist two opposite limits of the behaviour of this fifth force (though solutions anywhere in between could exist in regions and/or epochs of a real universe):

• When |f_R| ≪ |Φ|, we approximately recover the GR solution R = −8πGρ_m, which simplifies equation (5) to the usual Poisson equation:

  \begin{equation} \vec{\nabla }^2\Phi = 4\pi G\delta \rho _m, \end{equation}

  (6)

  where |$\delta \rho _m\equiv \rho _m-\bar{\rho }_m$| is the matter density perturbation. The effect of the fifth force is strongly suppressed in this regime, which is the consequence of what now is known as the chameleon mechanism (Khoury & Weltman 2004).
• In contrast, in the case of |f_R| ≥ |Φ|, we have |δR| ≪ δρ_m where |$\delta R\equiv R-\bar{R}$|⁠, leading to a modified Poisson equation:

  \begin{equation} \vec{\nabla }^2\Phi \approx 16\pi G\delta \rho _m/3, \end{equation}

  (7)

  which describes a 1/3 enhancement compared to the strength of Newtonian gravity. This enhancement is independent of the functional form of f(R), though the specific choice does affect the transition between this regime and the previous one. In the literature, these are often called the screened and unscreened regimes, respectively.

The chameleon mechanism earns its name because of the apparently environmental dependence of the fifth force; namely, that it is screened in regions where the Newtonian potential is deep enough. Intuitively, the scalar field f_R is the mediator of the fifth force, and has a mass due to its self-interaction (as described by its potential), so that the fifth force is of Yukawa type. The mass is environment dependent and becomes quite heavy in deep Newtonian potentials, causing the Yukawa force to decay exponentially with distance, which would make it difficult to detect experimentally the fifth force on Earth.

The Solar system is an example of a scale where screening could occur, making f(R) gravity viable (i.e. not yet ruled out by fifth force experiments). However, determining whether the screening has indeed taken place is much more tricky, because this also depends on the environment of the Solar system itself, such as the Galaxy and its underlying dark matter: even if the Solar system itself is incapable of producing a deep enough Newtonian potential to trigger screening, the Galaxy and its dark matter halo might just be able to do that.

In this paper, we shall not pursue this line of research. Instead, we will focus on cosmological scales and, in particular, regions of low matter density – the so-called cosmic voids. The idea is that even if a given f(R) model passes all local tests, it can still deviate substantially from GR in voids, where the chameleon screening does not happen. This opens a window of testing gravity independently.

2.1.2 The choice of f(R) model

The mere requirement of chameleon screening taking place in the Solar system for the model to be viable does not place a strong restriction on the functional form of f(R) which, as mentioned above, determines how the fifth force transitions from screened to unscreened regimes. Thus, in the literature many possibilities have been studied to different levels of detail.

However, there are good reasons for focusing on one particular example in a detailed study. First, at this point, neither the motivation for the functional form of f(R) nor the general understanding of the origin of the cosmic acceleration has matured to allow a connection to a more fundamental theory, and so any choice of f(R) is phenomenological. Indeed, f(R) models are known to face difficulties if inspected from a theoretical point of view (see e.g. Erickcek et al. 2013).

Secondly, although the exact details differ, many f(R) models share some common properties – chameleon screening in dense regions being one – resulting in qualitatively similar behaviour on cosmological scales. As a result, instead of asking ‘what do observations tell us about the form of f(R)?’, we can ask the question ‘to what extent are deviations from GR in the way that is prescribed by f(R) gravity allowed by observations?’ This question can be answered by choosing one particular model and testing it against as many observations as possible, as accurately as possible.

For these considerations, in this paper we choose to work with the model first proposed by Hu & Sawicki (2007, hereafter HS), which has three advantages. First, it is the most well-studied choice of f(R) gravity in the literature, so that this paper is building upon many existing tests of the model. Secondly, its particular f(R) form allows an efficient algorithm to be used to greatly speed up numerical simulations of the model (Bose et al. 2017), making its study easier, more efficient, and more accurate. Finally, not only is this model a representative example of f(R) gravity but it also – as we shall see shortly – has the flexibility to cover very different behaviours by varying its parameters. This last property means that studies of this particular model can inform us of general implications beyond the f(R) class of models.

In this paper, we will focus on the particular choice, n = 1, and consider two values of f_R0 – |f_R0| = 10⁻⁶, 10⁻⁵, which are referred to as F6 and F5 in the literature. These choices of f_R0 cover a broad range of model behaviours: F6 is the one that deviates least strongly from GR, while F5 deviates the most. Recent studies show that F5 could already be in tension with various cosmological observations (see e.g. Lombriser 2014 for a review of current constraints on f(R) gravity), but for completeness, we will still consider F5 in what follows.

2.2 Simulations and halo/galaxy catalogues

The simulations used in this paper (elephant simulations, or Extended LEnsing PHysics using ANalaytic ray Tracing) were performed using the ecosmog code (Li et al. 2012a), which is an augmented version of the popular simulation code ramses (Teyssier 2002), with added new routines to solve the scalar field and Einstein equations in modified gravity models. It inherits ramses's efficient mpi parallelization and adaptive-mesh-refinement (AMR) feature. The code starts with a uniform domain grid that covers a cubic simulation volume and the grid has |$N^{1/3}_{\rm dc}$| cells a side. The cells are refined (split into eight daughter cells) if the particle number in them exceeds some criterion (N_ref); in this way, the code hierarchically refines the domain grid to achieve high resolutions for solving the fifth force when the particle number density is high. The rest of the code algorithm is the same as in simulations for standard gravity, and interested readers are referred to Teyssier (2002, for more details).

The cosmological and numerical parameters are listed in Table 1, and the former are chosen to correspond to the best-fitting (Hinshaw et al. 2013) ΛCDM parameters. The simulations were initialized at z_ini = 49, with the initial conditions generated using the publicly available mpgrafic code (Prunet et al. 2008). As a control, for each realization of f(R) simulations, we also ran a ΛCDM simulation with the same parameters and specifications. Since the F5 and F6 models only deviate from ΛCDM at late times, at z_ini the effect of the fifth force is negligible, and so we let the simulations of all models start from exactly the same initial conditions. For ΛCDM, we had five independent realizations of boxes whose initial conditions differ only in their random phases, and we will use these to estimate errors in this paper. For f(R) gravity, only the first two realizations were run up to z = 0 and used for the analysis here.

2.2.1 Halo and galaxy catalogues

Dark matter haloes and the self-bound substructures associated with them are identified using the publicly available rockstar halo finder (Behroozi, Wechsler & Wu 2013).¹rockstar uses the six-dimensional phase-space information from the dark matter particles to identify haloes. Note that, in principle, the presence of the fifth force in f(R) gravity would require a modification to the unbinding procedure employed by rockstar (or indeed any halo-finding algorithm); however, Li & Zhao (2010) found the effect of this modification to be quite small and, as such, we use identical versions of rockstar for the GR and f(R) simulations. In this paper, we make use of only independent (‘main’) haloes, and not their substructures.

The parameters of the HOD are calibrated by requiring that the galaxy catalogue produced matches up with the galaxy distribution obtained from redshift surveys according to various metrics: most commonly, the number density of galaxies, n_g, and their projected two-point clustering. For our catalogues, the target (comoving) number density of galaxies were chosen to be appropriate for the boss cmass dr9 galaxy sample (Anderson et al. 2012): n_g = 3.8 × 10⁻⁴ h³ Mpc⁻³ and n_g = 3.2 × 10⁻⁴ h³ Mpc⁻³ at z = 0 and 0.5, respectively. Simply specifying the target number density is not enough to constrain the HOD, so we additionally require that the galaxy distribution exhibit the same two-point clustering across all gravity models at each redshift. Imposing this criterion in addition to the number density target substantially reduces the degree of degeneracy between different permutations of HOD parameter values. In total, we have constructed three types of HOD catalogues:

• cmass-default: in which the same set of HOD parameters are used for all gravity models;

• cmass-fixed-n_g: in which the parameters of the HOD are tuned (separately for each gravity model) to reproduce the same number density of galaxies only;

• cmass-fixed-n_g-ξ_g: in which the parameters of the HOD are tuned (separately for each gravity model) to reproduce the same number density of galaxies and projected clustering.

For catalogue (i), we use the same set of HOD parameters as in Manera et al. (2013), which were tuned to create a mock galaxy catalogue representative of the boss cmass dr9 observational sample at z = 0.5. The catalogues in (ii) were tuned separately at z = 0.5 and z = 0 to match the target number densities for that redshift. This was done by varying M_min while keeping the mass ratios, M₀/M_min and M₁/M_min, constant. The catalogues in (iii) were obtained by varying all five HOD parameters for each f(R) model, redshift, and realization so that it has the same galaxy number density and two-point correlation as the GR mock catalogue for the same redshift and realization. The tunning of the HOD parameters was achieved using a search with the Nelder–Mead simplex algorithm through the five-dimensional parameter space that minimized the rms difference between the f(R) and GR two-point correlation function in the distance range [2, 80] h⁻¹ Mpc (for more details, see Li & Shirasaki 2018). We tried different starting points for the search algorithm and all converged to the same value, suggesting that the choice of HOD parameters is reasonably unique. Table 2 summarizes the HOD parameter values used to construct the cmass-fixed-n_g-ξ_g catalogues for the first (Box-1) and second (Box-2) realizations. Note that for both z = 0 and z = 0.5, the GR parameter values are the same as those in Manera et al. (2013) (i.e. the same as catalogue i).

The difference in galaxy clustering induced by the different HODs is shown in Fig. 1, which compares the galaxy two-point correlation function, ξ_g, in GR, with the F5 model for catalogues (i)–(iii). The top panel shows that simply fixing the F5 HOD parameters to the same values as in GR (cmass-default), or fixing them to achieve the same number density (cmass-fixed-n_g) does not guarantee that the clustering is the same in the two models, and can in fact leave differences as large as 5 per cent across all separations. The bottom panel of Fig. 1 compares ξ_g in GR, F6, and F5 at z = 0 and z = 0.5, but for the cmass-fixed-n_g-ξ_g catalogue only; the agreement between the models now improves to better than 2 per cent for separations larger than 2 h⁻¹ Mpc. Fig. 2 displays the mean galaxy occupancy, 〈N(M)〉, as a function of halo mass, M, for the cmass-fixed-n_g-ξ_g catalogues in GR, F5, and F6. Finally, we note that catalogues (i) and (ii) are only used for illustration purposes in Fig. 1; the fiducial catalogue used in the rest of the analysis is catalogue (iii).

Figure 1.

Top panel: the two-point galaxy correlation function, ξ_g, in F5 for various z = 0 HOD catalogues. The solid black curves show the cmass HOD for GR. The other curves show the following F5 HOD catalogues: the same galaxy number density n_g and ξ_g as GR (cmass-fixed-n_g-ξ_g), only the same n_g as GR (cmass-fixed-n_g) and default GR HOD parameters (cmass-default). The secondary panel shows the ratio of the F5 ξ_g to the GR one. Bottom panel: ξ_g for the cmass-fixed-n_g-ξ_g GR and f(R) catalogues at z = 0 and z = 0.5. The secondary panel shows the ratio of the f(R) correlation function to the GR one. The grey shaded region show the 1σ uncertainty for the GR two-point correlation function.

Open in new tabDownload slide

Figure 2.

The mean number of galaxies, 〈N〉, as a function of halo mass, M, for our cmass-fixed-n_g-ξ_g catalogues. The different curves show the GR, F6, and F5 HOD models. The secondary panel shows the ratio of the f(R) HOD to the GR one.

Open in new tabDownload slide

2.2.2 Mock observations

We study the populations of voids at two redshifts, z = 0.0 and 0.5. The voids are identified using HOD galaxy catalogues that are built from the distribution of dark matter haloes in the simulation snapshots corresponding to the two redshift choices. For simplicity, when we apply the void finders to the distribution of galaxies in the periodic simulation box, which has a side length of 1024 h⁻¹ Mpc, we work in the distant observer approximation, which is when all the line of sights can be taken to be parallel, and we fix the line of sight to be along a preferential axis of the simulation box. This means that the two-dimensional (2D) underdensities, which are identified in the projected galaxy distribution, correspond to geometrical cylinders along the line of sight; in actual observations or in more realistic mock catalogues, the 2D underdensities would correspond to conical frustums (cones with their top sliced off). We work within the same distant observer approximation when calculating projected matter density profiles and the tangential shear signal. Our mock galaxy catalogues also neglect redshift space distortions, which, while not important for 2D voids, can affect the identification of 3D voids. A detailed study of the redshift-space distortions effect on void finding and void properties will be left for future work.

3 VOID FINDERS

Our study makes use of several void finders with the goal of determining which ones are best suited for probing chameleon-type modified gravity models. Broadly, we can split the void finders into two categories: the ones that identify 3D underdensities (i.e. voids) and the ones that identify 2D underdensities in the plane of the sky (2D spherical void finders, tunnels, and troughs). The underdensities identified by the latter methods are not formally called voids, since voids are 3D objects, but nonetheless, we can think of them as 2D voids that are very elongated along the line of sight. We identify the 2D underdensities by projecting the full simulation box along one of the principal axes of the simulations. The box has a side length of 1024 h⁻¹ Mpc, which is approximately the comoving distance between redshift 0.3 and 0.7 – the latter is the redshift range we use to make predictions for euclid- and lsst-like surveys (see Section 4.5). The structures identified by each method are illustrated in Fig. 3. The rest of this section provides a short description of each of the six methods used to classify 2D and 3D underdensities.

Figure 3.

Illustration of the structures identified by the six void-finding methods employed in this paper. Each circle corresponds to an underdensity of radius equal to the one shown in the plots. The left column plots the svf (top left), wvf (centre left), and zobov (bottom left) 3D voids in a 50 h⁻¹ Mpc slice, with the background image showing the density in that slice. The right column plots the 2D svf_2d voids (top right), tunnels (centre right), and troughs (bottom right), with the background image showing the projected density of the full box (which has a 1024 h⁻¹ Mpc side length) along the line of sight. Note the different scales for the left and right columns.

Open in new tabDownload slide

3.1 3D spherical underdensity void finder

The 3D spherical underdensity void finder (from hereon svf) used in this work is a substantially modified version of the algorithm presented in Padilla et al. (2005). It searches for spherical regions that satisfy a specified density criterion in a simulation.

To find the void centres, a rectangular grid is constructed over the simulation volume. The number of galaxies in each cell of the grid is counted, and the centres of the cells that are empty of galaxies are considered to be prospective void centres. The integrated density profile about each centre is then calculated. The radius at which the integrated density is equal to 20 per cent of the mean number density of galaxies is considered to be the radius of a void. Only the largest sphere satisfying this condition about any centre is kept in the void catalogue. This density threshold is commonly adopted in the literature of void studies, motivated by a calculation presented in Blumenthal et al. (1992), who used linear theory to show that the voids we observe at the present time correspond to mass underdensities of about 20 per cent of the cosmic mean. After this step, one might end up with a large number of voids with similar centres, sharing a fraction of their volume. If the distance between the centres of any two voids is less than 80 per cent of the sum of their radii, we keep only the largest of the two. This naturally favours a single, large void over many smaller voids of comparable sizes and similar centres. A conservative degree of overlap between adjacent voids is still required, though, in order to trace the low density regions of the cosmic web to a good extent.

There are two free parameters in the 3D spherical underdensity void finder, the density criterion and the overlapping threshold for excluding voids. As mentioned above, the former is chosen as 20 per cent of the mean number density of galaxies for physical considerations. For the latter, Cai et al. (2015) have tested several different choices and found that the S/N is not sensitively dependent on the exact value, and so we choose 80 per cent as our default.

3.2 3D Watershed void finder

The watershed void finder (wvf; Platen et al. 2007) associates the voids to the watershed basins of the large-scale density field without imposing a priori constraints on the size, shape, and the mean underdensity of the objects it identifies. The method starts by constructing a galaxy density field using the Delaunay Tessellation Field Estimator (Schaap & van de Weygaert 2000; Cautun & van de Weygaert 2011), which uses a Delaunay triangulation with the galaxies at its vertices to extrapolate a volume-filling density field. The resulting density is defined on a 1024³ regular grid with a grid cell size of 1 h⁻¹ Mpc. The density is then smoothed with a Gaussian filter of 2 h⁻¹ Mpc radius to reduce small-scale structures inside and at the boundaries of voids, which could potentially give rise to artificial voids. The 2 h⁻¹ Mpc filter corresponds to the typical width of the filaments and walls forming the void boundaries (e.g. Cautun et al. 2013, 2014). The smoothed density field is then segmented into watershed basins. This process is equivalent to following the path of a rain drop along a landscape: each volume element, in our case the voxel of a regular grid, is connected to the neighbour with the lowest density, with the same process repeated for each neighbour until a minimum of the density field is reached. Finally, a watershed basin is composed of all the voxels whose path ends at the same density minimum.

The resulting wvf void catalogue is characterized in terms of the void centres, which are chosen as the volume-weighted barycentre of all the voxels associated with each void, and the void sizes. Since voids have irregular shapes, the latter is given in terms of the effective void radius, |$R_{\rm eff} =(\frac{3}{4\pi }V)^{1/3}$|⁠, which is the radius of a sphere with the same volume as the void volume.

The wvf has a single free parameter – the Gaussian smoothing size. As mentioned above, our choice of 2 h⁻¹ Mpc corresponds roughly to the size of filaments and helps to reduce artificial voids. The smoothing size is much smaller than the radius of most of the wvf voids considered in this paper, and we do not expect it to strongly affect the properties of the wvf voids.

3.3 3D zobov void finder

We also run a slightly modified version of the zobov algorithm (Neyrinck 2008) on our HOD galaxy mocks. zobov uses Voronoi tessellations to estimate the galaxy density field at each galaxy position and then identifies the density minima by comparing each Voronoi cell with its neighbours. Starting from the density minima, neighbouring Voronoi cells of increasing densities are grouped together to form ‘zones’. The growth of ‘zones’ stops when the density of the next neighbouring Voronoi cell decreases. The original algorithm goes one step further to group neighbouring ‘zones’ together if their boundary is below a density threshold, i.e. as this threshold density increases, more and more ‘zones’ are grouped together to form larger and larger voids. This step can lead to unwanted effects; Cai et al. (2017) showed that when applying ‘zone’ merging some of the largest zobov voids found in their cmass mocks do not seem to correspond to true matter underdensities. To avoid such spurious effects, we treat every ‘zone’ as a void in this work, i.e. we do not merge zones. By definition, our zones (voids) do not overlap with each other.

The void volume, V, is defined as the sum of the Voronoi cell volumes that are part of that given ‘zone’. The effective radius of the void is then defined as |$R_{\rm eff} =(\frac{3}{4\pi }V)^{1/3}$|⁠. Voronoi cells belonging to each zone are weighted by their volumes to define the void centre.

Our implementation of the zobov void finder has no free parameters given a tracer catalogue.

3.4 2D Spherical underdensity void finder (svf_2d)

2D spherical voids are obtained using a slightly modified version of the svf algorithm. To find the void centres, a rectangular grid is constructed over the projected distribution of HOD galaxies along one of the axes of the simulations, and the number of galaxies in each grid cell is counted. The centres of empty grid cells are considered as prospective void centres, and circles are grown from those centres until the integrated number density of galaxies at the circle radius is equal to 40 per cent of the mean density. After this step, if two voids overlap more than 80 per cent of the sum of their radii, only the largest void is kept in the catalogue.

The svf_2d, as the 3D one, has three parameters: the density criterion for defining a void, the criterion for removing overlapping voids, and the length of projected redshift range. For the former, we have tested density criteria equal to 0.2, 0.3, 0.4, 0.5, and 0.8 of the mean projected galaxy number density, respectively, and found that choosing the 0.4 value results in the strongest detection of weak lensing by 2D underdensities (when considering only sample variance uncertainties). Similar to the 3D case, if the centres of two neighbouring voids are closer than 80 per cent of the sum of their radii, we remove the smaller one. Here, we identify voids by projecting in the distant observer approximation the entire simulation box (1024 h⁻¹ Mpc in length) along one of its preferential axes.

3.5 2D tunnels

The tunnels correspond to elongated line-of-sight regions that intersect one or more voids without passing through overdense regions (Cautun et al., in preparation). Using galaxies as tracers of the matter distribution, the tunnels are identified as circles in the plane-of-the-sky that are devoid of galaxies. The typical size of tunnels depends on the number density of tracer galaxies and on the line-of-sight depth used to identify them; a higher tracer density or a larger line-of-sight depth results in smaller tunnels. In the distant observer approximation, which we use in this work, the tunnels consist of line-of-sight cylinders that are empty of galaxies.³

To identify the tunnels, we start by projecting the HOD galaxy catalogue along one of the simulation axes to obtain a 2D distribution of galaxies. To identify the largest circles empty of galaxies, we build a Delaunay tessellation with the galaxies at its vertices. By definition, the circumcircle of every Delaunay triangle is empty of galaxies, with the closest galaxies being the three that give the triangle vertices and that are found exactly on the circumcircle. The tunnels consist of the circumcircles whose centres are not inside a larger circumcircle. We also discard candidates for which the Delaunay triangle has an area of 0.2 or less than that of its circumcircle; such cases correspond to triangles that either have a side much shorter than the other two or that have one very large angle. The tunnel centre and radius are given by the centre and radius of its corresponding circumcircle.

Since galaxies are biased tracers of dark matter, tunnels do not necessarily correspond to empty regions of dark matter. Nevertheless, the tunnel radius is correlated with their projected matter density, as shown in Appendix A. On average, large tunnels correspond to underdense regions while small ones correspond to overdense ones, with the transition taking place at a tunnel radius of 1 h⁻¹ Mpc, which corresponds to 0.4 times the mean projected galaxy separation. Here, we only consider tunnels larger than this transition radius, since we are interested in the modified gravity signature of underdense regions.

The tunnel catalogue depends on two free parameters, the transition radius between underdense and overdense tunnels, and the line-of-sight projection length. The former is determined by analysing the enclosed projected matter density inside tunnels. The latter, which affects the size distribution of svf_2d voids too, should be selected according to the details of the observational survey we want to match. For this study, we project the entire simulation box (1024 h⁻¹ Mpc in length) using the distant observer approximation.

3.6 2D troughs

The troughs (Gruen et al. 2016) are somewhat similar to tunnels in that they also correspond to elongated line-of-sight regions that are underdense. In contrast to tunnels, all troughs have the same radius and are given by randomly positioned circles in the plane of the sky that contain very low galaxy counts.

The troughs are identified using the same projected HOD galaxy catalogues as the tunnels, which represent our simulated plane of the sky in the distant observer approximation. We choose to study troughs of 2 h⁻¹ Mpc in radius, which is similar to the typical radius of tunnels. Moreover, troughs of this size contain a similar galaxy count at the 5 arcsec troughs studied by Gruen et al. (2016) and Barreira et al. (2017). The trough identification starts by randomly positioning 10⁶ circles of 2 h⁻¹ Mpc in radius on our simulated plane of the sky. The troughs are given by the circles that contain two or fewer galaxies inside them; these correspond to 23 and 30 per cent of population at z = 0 and z = 0.5, respectively. This selection is similar to the one used by Gruen et al. (2016), and it is also the one that gives a good compromise between selecting very underdense regions and covering a large fraction of the available simulation. See Appendix A for a more detailed discussion.

3.7 Comparison of the different void finders

Fig. 3 presents a visual comparison of the underdensities identified by the six void finders. For the 3D voids, the figure shows the voids whose centres lie in 500 × 500 × 50 (h⁻¹ Mpc)³ region of the simulation. The background image shows the density, ρ, of that region expressed in terms of the mean background density, |$\overline{\rho }$|⁠. While some voids seem to have overdense regions inside them, in most cases it is either a projection effect or it is due to representing the wvf and zobov voids as circles when they typically have non-spherical shapes. A closer inspection of the positions of 3D voids reveals that, in a large number of cases, we can find a match between individual objects identified by different void finders. This is especially striking for the wvf and zobov methods, since both are based on the watershed transform of the galaxy density field. Nonetheless, many of the matched voids have different centres and radii, and thus no two methods result in similar void catalogues (see Colberg et al. 2008 for a comparison of more void finders).

The 2D underdensities generally have much smaller radii than their 3D counterparts, so the right column of Fig. 3 shows 2D voids in a smaller region of size 25 × 25 (h⁻¹ Mpc)² in the xy plane of the whole simulation box. The 2D voids were identified by projecting the full simulation box (1024 h⁻¹ Mpc side length) along the z-direction in the distant observer approximation, and the background colours show the projected matter density field (in unit of mean projected density) along the z-direction for the full simulation box. The three right-hand panels show that the 2D voids, which were selected as underdensities in the projected galaxy distribution, correspond to underdensities in the projected matter density field too. Most matter underdense regions are identified as 2D voids, but the centre and size of those 2D voids show a large variation between the three types of objects we study: svf_2d voids, tunnels, and troughs. Both svf_2d voids and tunnels show relatively little overlap, but on average the tunnels are smaller and cover a larger fraction of the plane of the sky than svf_2d voids. The troughs are very clustered in underdense regions and show a large degree of overlap; in fact, for readability purposes, the bottom right-hand panel of Fig. 3 shows only half of the troughs we identified, with the half not shown being right on top of the ones shown in the figure.

4 RESULTS

In this section, we compare the distribution of void properties, e.g. void abundance, galaxy and matter density profiles, and lensing signal, between the standard ΛCDM cosmology and f(R) models. We perform this analysis using the cmass-fixed-n_g-ξ_g HOD catalogues, which means that both the GR and f(R) mocks have the same number density of tracers and the same real-space two-point correlation function of these tracer galaxies. We perform this comparison for voids defined using six different methods: three that identify 3D underdensities (svf, wvf, and zobov) and three that identify 2D underdensities using the projected distribution of galaxies (tunnels, troughs, and 2D wvf).

We characterize voids in terms of their radius, which we denote with R_eff and R_p eff (i.e. projected radius) for the 3D and 2D structures, respectively. This nomenclature is motivated by two of the void finders, wvf and zobov, that identify irregularly shaped voids and thus these voids are characterized in terms of an effective radius, which is the radius of a sphere with the same volume as that of the void. For consistency, we use the same R_eff notation also for the radius of the other voids, even though those are by definition spherical/circular objects.

We present stacked void profiles that are averaged over voids of all sizes, unless we specify otherwise. We have checked that similar differences between GR and f(R) are present when stacking voids in a narrower range of sizes, albeit with a lower S/N due to reduced sample size and thus poorer statistics. When averaging, each void is given the weight, |$w= R_{\rm eff} ^2$|⁠, such that larger voids are weighed more compared to the more numerous small voids. This weight is motivated by measurements of the tangential shear profile of voids in observations. The number of source galaxies inside a void scales approximately with |$R_{\rm eff} ^2$| and thus large voids have a larger contribution to the stacked tangential shear measurements. We note that when simulating other observables, different weights may be appropriate; for example, in the case of the 3D void galaxy density profile, the contribution of each void to the stacked signal scales with |$R_{\rm eff} ^3$|⁠. Nonetheless, for simplicity, we keep the same |$w= R_{\rm eff} ^2$| weight for all our stacked profiles.

We present mean profiles by averaging over two realizations of a cosmological volume of 1024 h⁻¹ Mpc in length, Box 1 and 2, which were simulated using both GR and f(R) gravity models. We estimate uncertainties using five realizations of the GR box. For each realization, we split the volume into 4³ non-overlapping regions and perform 100 bootstrap samples over these regions. The procedure leads to 5 × 100 samples that we use to compute the correlation matrix and estimate errors. For tunnels and troughs, we split the mock plane of the sky into 8² non-overlapping regions, after which we perform exactly the same analysis. Such a procedure accounts for correlations between neighbouring voids that are especially important in the case of troughs since many troughs overlap. In practice, the procedure is implemented by first computing the stacked profile and the total weight of the voids in each of the regions defined above. The total weight of all voids inside a region is assigned as the weight of that region. Then, all the regions selected for a bootstrap sample are combined according to their weights.

4.1 Void abundances

We start by comparing the distribution of void sizes between GR and f(R) models. Fig. 4 shows these distributions for all the void finders employed here with the exception of troughs that are defined to have a constant 2 h⁻¹ Mpc radius. The distribution of void sizes varies between finders (note that the horizontal axes have different ranges in the different panels of Fig. 4), but, when considering the same identification method, there is no statistically relevant difference between GR and f(R) gravity. This result agrees with that in Cai et al. (2015), where it was found that while voids identified in the matter density field are larger in f(R) models compared to the GR case, these differences largely disappear when identifying voids using the halo distribution. Falck et al. (2018) obtained a similar result for the DGP modified gravity model: when identifying voids using haloes, there is no difference in void abundance between ΛCDM and DGP models. In contrast, Zivick et al. (2015) found that f(R) models boost the number of large voids; the discrepancy could be due to their neglect of halo and galaxy bias since they identified voids using a subsampled distribution of dark matter particles (e.g. see Pollina et al. 2016).

Figure 4.

Comparison of the void abundance, i.e. number density of voids as a function of their radius, in GR and f(R) models using HOD galaxy catalogues that were tuned to have the same number density and two-point correlation function across all models. The figure shows the abundance of 3D voids identified using svf (top left-hand panel), wvf (centre left-hand panel), and zobov (bottom left-hand panel); and that of 2D svf_2d voids (top right-hand panel), and tunnels (centre right-hand panel). All the 2D troughs have by definition the same 2 h⁻¹ Mpc radius. The 2D voids were obtained by projecting the entire simulation box (1024 h⁻¹ Mpc side length) along one of its preferential axes in the distant observer approximation. For each panel, the various curves show the results for GR and for two f(R) models at z = 0 (solid curves) and z = 0.5 (dashed curves). For clarity, only the GR z = 0.5 are shown in each of the main panels. The secondary panels shows the ratio of the f(R) results to the GR one for both z = 0 and z = 0.5. The shaded region shows the 1σ uncertainty interval computed using multiple GR realizations.

Open in new tabDownload slide

By analysing our various HOD catalogues (see Section 2.2.1), we find that the void abundance is most sensitive to the number density of tracer galaxies, n_g. Once the f(R) galaxy catalogues have the same n_g as the GR ones, there is little difference in void abundance. Further matching also the two point correlation function of the f(R) and GR catalogues does not lead to a significant change. We have checked explicitly that (not shown here) the same is true if dark matter haloes with fixed number density are used as tracers to find voids. The differences in void abundance between z = 0 and z = 0.5 (the solid versus dashed curves from Fig. 4) could also be due to differences in galaxy number density, which are n_g = 3.8 × 10⁻⁴ and 3.2 × 10⁻⁴ h³ Mpc⁻³ at z = 0 and 0.5, respectively. A lower tracer number density results in systematically larger voids and fewer small voids due to them merging to form larger ones.

4.2 Void galaxy number density profiles

Fig. 5 presents the radial galaxy number density profiles of voids. All the methods find that, on average, the void interiors are devoid of galaxies, but the average galaxy number varies between methods. For svf, n_g shows a large increase at the void radius since svf voids are identified as the largest sphere that has mean enclosed density |$\le 0.2\overline{n_{g}}$|⁠, with |$\overline{n_{g}}$| denoting the mean number density of galaxies. The watershed void finders, both wvf and zobov, have smoothly increasing n_g(r) profiles due to their non-spherical shape, with their overdense ridges smeared by the spherical stacking procedure (Cautun, Cai & Frenk 2016). The profiles of 2D voids are expressed in terms of the projected mean galaxy number, Σ_g, and its mean background value, |$\overline{\Sigma _{\rm g}}$|⁠. Similar to the 3D svf voids, the svf_2d voids are very underdense inside the void and show a prominent galaxy density enhancement at their edge. By definition, the tunnels are the largest circles devoid of galaxies; this explains why Σ_g = 0 for r < R_p eff is followed by a sharp peak at r = R_p eff. The troughs have |$\Sigma _{g}<\overline{\Sigma _{g}}$| not only inside R_p eff but also outside their radius. This is because troughs with our selection criteria are likely to reside in large underdense regions⁴ (see Fig. 3). With the exception of troughs, the n_g(r) and Σ_g profiles show very little dependence on redshift.

Figure 5.

Comparison of the radial galaxy density profile of voids between GR and f(R) models. The left column shows the 3D galaxy number density profile, |$n_{\rm g}/\overline{n_{\rm g}}$|⁠, of 3D voids. The right column shows the projected surface density of galaxies, |$\Sigma _{\rm g}/\overline{\Sigma _{\rm g}}$|⁠, of 2D line-of-sight underdensities. The symbols and curves are the same as in Fig. 4.

Open in new tabDownload slide

Fig. 5 shows that there is no statistical difference in the galaxy number density profiles between GR and f(R) models when the latter catalogues are matched to have the same number density of tracers and the same galaxy correlation function. The only exception is inside svf_2d voids and troughs, where the differences are likely a limitation of our mocks rather than a tell-tale signature of modified gravity. This is supported by the non-monotonic behaviour with redshift and with the parameter that determines the strength of the modifications to gravity: for example, the F5 model, which has stronger deviations from GR than F6, shows a large systematic difference at z = 0.5 but no difference at z = 0, while in theory we would expect bigger differences at later times. For F6, it is especially suspicious that the z = 0 difference in the profile with respect to GR is bigger than the difference of the F5 model. The observed differences in the galaxy density profiles within the trough radius, which is 2 h⁻¹ Mpc, could be due to the fact that our HOD catalogues were tuned to match the galaxy correlation functions among the different models only above separations of 2 h⁻¹ Mpc. In particular, we note the same qualitative (and even quantitative) behaviour of the model differences in the case of 2D spherical voids, suggesting that this is not specific to troughs. The error bars in the plot appear smaller in the case of troughs, but these are simply taken as the square root of the diagonal elements of the covariance matrix for indication, and the actual errors in different bins are indeed very strongly correlated. Therefore, we refrain from interpreting the differences shown in the galaxy density profiles of troughs and svf_2d voids as a signature of deviation from GR, and will leave a more detailed investigation to a future work, hopefully using higher resolution simulations and mock catalogues.

4.3 Void matter density profiles

The resulting matter density profiles of voids are shown in Fig. 6. These are obtained by stacking all the voids and thus are average profiles. All methods identify regions that are underdense in the inner parts, i.e. r ≲ R_eff; and that show a modest overdensity close to their edge, i.e. r ≃ R_eff. The only exception are troughs that are underdense even beyond r > R_eff because troughs represent the most underdense regions and moreover, because their boundaries are not defined as a galaxy number overdensity. Compared to the galaxy number density profiles, the void matter profiles are both less underdense in their interiors and less overdense at their boundaries. Moreover, the matter profiles show a considerable dependence with redshift, with voids identified at higher redshift being less overdense.

Figure 6.

Comparison of the void matter density profile between GR and f(R). The left column shows the 3D matter density profile, |$\rho /\overline{\rho }_{bg}$|⁠, of 3D voids. The right column shows the projected matter surface density, |$\Sigma /\overline{\Sigma }_{bg}$|⁠, of 2D line-of-sight underdensities. The symbols and curves are the same as in Fig. 4.

Open in new tabDownload slide

We find that the interiors of both 2D and 3D voids are more underdense in f(R) gravity than in GR. The effect is the strongest for F5 and also decreases at higher redshift. This is in line with theoretical expectations, since the modifications to gravity become stronger at later times and are larger in F5 than in F6. Due to mass conservation, since f(R) voids have emptier interiors they should also be more overdense around the void edge. This effect is present for all void types, except troughs, and it is especially prominent for svf voids, both 2D and 3D ones, and tunnels. The troughs are different, since they are likely embedded in underdense regions much larger than their size, and, even by going to 2R_eff, we still have not reached the overdense ridge surrounding these underdense regions.

We studied the void density profiles of several modified gravity HOD catalogues (see Section 2.2.1): cmass-default, which uses the same HOD parameter values as GR, cmass-fixed-n_g, which changes the three HOD mass parameters proportionally to match the galaxy number density of GR, and cmass-fixed-n_g-ξ_g (the one shown in Fig. 6), which matches both the number density and the two-point clustering of GR. For all these catalogues, the void density profiles show the same difference between f(R) and GR, suggesting that the difference between the models is robust to small changes in how galaxies populate dark matter haloes.

4.4 Void tangential shear profiles

Weak gravitational lensing represents the most promising way of measuring the matter distribution in and around voids, and could potentially be used as a probe of modifications to gravity (Barreira et al. 2015; Cai et al. 2015). Current surveys have already measured the weak lensing imprint of voids (e.g. Gruen et al. 2016; Sanchez et al. 2017) and future surveys, thanks to increases in both sky coverage and image quality, would improve greatly the precision of such measurements. This motivates us to compare the weak lensing signal of the various void finders and predict which method shows the largest potential for discriminating between GR and f(R) models.

According to equation (17), up to a normalization constant, the tangential shear is determined by the differential surface mass density, ΔΣ(r), which, in turn, is determined by the projected surface mass density at r, Σ(r), whose calculation is described in Section 4.3. In the case of 2D voids, Σ(r) is shown in the right column of Fig. 6. Similarly, we compute Σ(r) for the case of 3D voids by projecting the voids on the simulated plane of the sky using the distant observer approximation.

In the case of 3D voids, the small amplitude of the tangential shear signal introduces an additional challenge. This is illustrated in Fig. 7, where we show ΔΣ(r) for wvf voids in the GR model. The differential surface density was calculated as described above, by projecting the void and matter distribution along the line of sight. The solid curve shows one mock observation (i.e. one simulation box and one line-of-sight choice), while the uncertainty region corresponds to the 1σ sample variance calculated using multiple realizations and line of sights. The 1σ uncertainty is around 10 per cent of the signal strength at r ∼ R_eff and is much larger than the typical difference in the void density profile between f(R) and GR models, which is only a few per cent (see Fig. 6). This means that to systematically characterize the lensing differences between f(R) and GR, we would need a large number of mocks.

Figure 7.

Comparison of the differential surface mass density, ΔΣ, of wvf voids measured directly from the projected particle distribution and the one calculated from the 3D density profile using equation (20). Both methods were applied to the same physical volume. The uncertainty associated with the red curve is roughly the width of the curve itself and thus is barely visible. Calculating ΔΣ using the void 3D density profile gives a better estimate of the mean signal, but it underestimates the errors by a factor of ∼7, since this method does not account for a major source of error: the variation in the projected mass distribution of individual voids with the line of sight along which the void is observed. The blue-shaded region accounts for this effect and is the observationally relevant uncertainty.

Open in new tabDownload slide

A second curve in Fig. 7 shows the outcome of calculating the differential surface mass density of wvf voids using equation (20). As expected, we find good agreement between the mean values of the two calculations. The small discrepancy at r > 1.5R_eff is a combination of correlated errors and a slight overestimation of ΔΣ(r) due to the fact that we limit the integration in equation (20) to three times the void radius. The sample variance 1σ uncertainty region, which was found using multiple realizations of the GR box combined with bootstrap sampling, is much smaller than the uncertainty resulting from projecting the voids and matter distribution along the line of sight. This suggests that integrating the 3D matter density profile is a more computationally efficient method of calculating the mean ΔΣ(r) value of 3D voids. However, the same calculation underestimates the size of the sample variance error by a factor of ∼7. This discrepancy is due to equation (20) neglecting one major source of scatter. The ΔΣ(r) uncertainty is a combination of two effects. First, it is affected by the void-to-void variation in their radial mass distribution. Secondly, since voids are highly non-spherical, the projected mass distribution around each void shows plenty of variation depending on the viewing direction. The 3D density profile of each void corresponds to an average over all possible line of sights and thus does not include this latter source of scatter. Therefore, the observationally relevant uncertainty is the one computed using the projected particle distribution.

To summarize, for 3D voids we follow a hybrid approach for calculating ΔΣ. The mean ΔΣ for both GR and f(R) models was computed using equation (20), that is by integrating the 3D density profiles along the line of sight. This means that we can measure to a high accuracy systematic differences between GR and f(R) models. To compute the observationally relevant GR sample variance, we used the projected particle distribution around each void. Similarly to previous error estimates, we compute the covariance matrix using 500 bootstrap samples constructed from five simulation boxes (for details, see the fourth paragraph in Section 4).

Fig. 8 shows the differential surface mass density, ΔΣ, for the six voids studied here. In all cases, we find that ΔΣ is negative at least up to, r ≲ 1.5R_eff, indicating that voids, due to their underdense interiors, produce divergent lensing, which is similar to a concave lens. In contrast, high-density regions (e.g. clusters) give rise to convergent lensing, which is similar to a convex lens. For most void types, the diverging weak lensing signal peaks at the void radius, r = R_eff, with troughs being the exception for which the signal peaks at r ≃ 1.2R_eff. Of the 3D underdensities, svf voids produce the strongest tangential shear, which is |${\sim }30\, {\rm per \, cent}$| larger than the signal of the other two 3D voids. The lensing signature of wvf voids, and probably that of zobov ones, can be increased by a factor of ∼2 by stacking with respect to the boundary of these non-spherical objects (Cautun et al. 2016), which would result in a stronger lensing signal than the svf one. The 2D underdensities have an even stronger weak lensing imprint than the 3D ones, with troughs producing an approximatelyfivetimes larger signal than svf voids. The svf_2d voids and tunnels have an even larger tangential shear, roughly 15 times larger than that of svf voids. This is due to the fact that 2D voids are much smaller than 3D ones and therefore they are measuring the matter fluctuations at much smaller scales.

Figure 8.

Comparison of the differential surface mass density, ΔΣ, profile of voids in GR and f(R) models. These were calculated by projecting the mass distribution along one of the axes of the simulation box. The grey-shaded regions show the 1σ uncertainties corresponding to the sample variance of the GR signal. For 2D voids, the errors correspond to the sample variance of a (1.024 h⁻¹ Gpc)³ volume, which is the volume of each of our five simulation boxes. For 3D voids, the errors are very large so we show the sample variance for an euclid-like survey with a 10 (h⁻¹ Gpc)³ volume, which is 9.3 times larger than the volume of each of our simulation boxes. The symbols and curves are the same as in Fig. 4.

Open in new tabDownload slide

Fig. 8 also compares the ΔΣ(r) profiles in GR and f(R) gravity. Voids in modified gravity models show a stronger lensing signal than in GR, with the enhancement being largest for the F5 model at late times. To assess the significance to which these differences can be measured, Fig. 8 shows as a grey-shaded region the 1σ sample variance for GR. For 3D voids, the uncertainty corresponding to the volume of the simulation box is very large, so we present a rescaled error for an euclid-like survey with a 10 (h⁻¹ Gpc)³ volume (9.3 times larger than our simulation box), which practically means reducing the uncertainty by a factor of |$\sqrt{9.3}$|⁠. Among the 3D voids, svf voids show the largest difference compared to the fiducial GR model, both in absolute terms as well as when compared to the error bars, but note that the |$\Delta \Sigma _{f(R){}}-\Delta \Sigma _{\small {GR}}$| difference is very similar to the cosmic variance. In contrast, for 2D voids, the signature of modified gravity models is significantly larger than the associated uncertainties in the lensing signal, which makes these objects ideal for testing f(R) models.

4.5 Predictions for euclid and lsst

Here, we investigate the extent to which void lensing in future surveys can be used to constrain modified gravity theories. In particular, we focus our attention to the euclid (Laureijs et al. 2011) and the lsst (LSST Science Collaboration et al. 2009) lensing surveys that cover 20 000 and 18 000 deg², respectively. For both surveys, we assume that their sky coverage area overlaps with spectroscopic surveys that have a galaxy number density in the redshift range z = 0.3–0.7 at least as high as that of the sdss cmass sample. In many cases, there will be overlapping spectroscopic surveys with higher galaxy number densities, in which case our analysis quantifies the modified gravity constraints that can be inferred when using only the brightest galaxies. Using the full spectroscopic survey would probably result in even tighter constraints, but our simulations lack the resolution to make predictions for such observations.

Here, we adopt z_s; med = 0.8 and 1.2, which corresponds to the median source redshift for the euclid and the lsst surveys, respectively, to obtain Σ_c; eff = 6770 and |$3960\,{h}\rm {M}_{{\odot }} \, \rm {pc}^{-2}$|⁠. For the euclid survey, we predict tangential shear values at the position of the dip of γ_t ≃ −1 × 10⁻⁴ for 3D voids, γ_t = −5 × 10⁻⁴ for troughs and γ_t = −2 × 10⁻³ for svf_2d voids and tunnels. For the lsst survey, we predict a lensing signal that is a factor of 1.7 times larger than for euclid. Thus, depending on the method used to identify underdense regions, the weak lensing signal can vary by a factor of 20, being lowest for 3D underdensities and highest for 2D underdensities, with svf_2d voids and tunnels having the highest lensing imprint.

The tangential shear measurements are affected by three important sources of uncertainty: void sample variance, the covariance of uncorrelated large-scale structure along the path of the light rays, and shape noise (e.g. see Krause et al. 2013). The first two error sources can be obtained by calculating the void tangential shear profile due to the mass distribution between the source plane and the observer. For this, we construct a mock light cone for each GR realization. First, the mass distribution in the redshift range z = 0.3 to z = 0.7 is given by the z = 0.5 snapshot of the respective GR realization. To account for uncorrelated large-scale structure, the mass distribution for z < 0.3 and for z > 0.7 is taken from the z = 0.5 snapshot of the other GR realizations. For simplicity, our light cone mocks use the mass distribution at z = 0.5, which neglects the time evolution of the clustering, and, secondly, we use a cylindrical geometry while in practice observations have a conical geometry. We calculate the tangential shear of each individual void by applying equation (17) to thin slices along the line of sight and summing the contribution of all these slices. We do so for all the 5 GR realizations.

The remaining source of uncertainty, shape noise, comes from the intrinsic ellipticity distribution, characterized by its variance, σ_ε, of the source galaxies used to measure γ_t. We measure shape noise by generating a random distribution of source galaxies in our simulated plane of the sky (i.e. projected simulation box), with each source having a randomly assigned and randomly oriented ellipticity, with the ellipticity variance being given by σ_ε (this is similar to the procedure described in Sanchez et al. 2017, but applied to mock catalogues and not to the data). For each void in the catalogue, we calculate the mean source galaxy ellipticity for the same radial bins used to estimate the weak lensing signal. For this calculation, we adopt and intrinsic source ellipticity, σ_ε = 0.22, and a number density of sources, n_sources = 30 and |$40\, \rm {arcmin}^{-2}$| for euclid and lsst, respectively. Then, we add the shape noise to the tangential shear signal for each of the voids found in the 5 GR realizations. To compute the covariance matrix, we split each GR void catalogue into 64 subregions and generate 100 bootstrap realizations; the resulting N = 500 estimates are used to calculate the total covariance matrix. The resulting covariance matrix, which corresponds to a survey with the same volume as each of our simulation boxes, is rescaled by multiplying with the factor V_sim/V_survey, where V_sim and V_survey are the comoving volumes of the simulation box and the survey, respectively.

We compute the cumulative S/N for all the six void catalogues for both the F5 and F6 models. The results are shown in Fig. 9. For all methods, the cumulative S/N increases up to the void radius, after which it stays relatively constant, which suggests that most of the power for distinguishing gravity models comes from the region r ≲ R_eff. The 2D voids show the largest S/N for both f(R) models, with tunnels being the most promising method. The S/N is larger for lsst than for euclid, since the source galaxies of the former survey are at higher redshift and there are 25 per cent more of them. In the case of euclid, tunnels have a maximum S/N of 50 and 7 for, respectively, the F5 and F6 models. For lsst, the tunnels’ S/N of the same two models increases to 80 and 11, respectively.

Figure 9.

The cumulative (from small to large radius) S/N of the differences in tangential shear between f(R) and GR. The top and bottom panels show the S/N of the F5 and F6 models, respectively. Each colour corresponds to one of the six void types studied in this paper. The predictions for an euclid- and a lsst-like lensing surveys are shown as solid and dashed lines, respectively.

Open in new tabDownload slide

For the F5 case, the lensing signal of 3D voids peaks at an S/N ∼ 1–2. This value is lower than previous studies, with Cai et al. (2015) predicting that svf voids in a 1 (h⁻¹ Gpc)³ volume can distinguish F5 with S/N ∼ 7. Extrapolating this result to the euclid volume would result in a |$\sqrt{10}$| larger value, i.e. S/N ∼ 22, while in this work, we only find an S/N ∼ 2. The discrepancy is due to several differences between our analysis and the Cai et al. one. We estimate the S/N using the z = 0.5 snapshots, while Cai et al. used the z = 0 matter distribution; as seen in Fig. 8, the difference in void tangential shear between f(R) and GR is smaller at higher redshift. Also, we include two additional error sources, shape noise, and the contribution of uncorrelated line-of-sight large-scale structure, which are similar in magnitude to the sample covariance of voids.

In Fig. 10, we study which of two possible observational strategies, which is surveying a larger volume versus having a higher density of source galaxies, optimizes the gain of using void lensing to distinguish between the f(R) and GR models. We analyse the case of F6 and we limit the analysis to the lsst survey, which in Fig. 9 gives the largest S/N. For each void finder, Fig. 10 shows three sets of S/N curves: (i) lsst (solid lines), (ii) four times the lsst volume, but equal number of lensing source galaxies as lsst (dotted lines), and (iii) lsst volume with four times the number of lsst lensing source galaxies. Cases (ii) and (iii) correspond to decreasing separately the covariance matrix contribution from, respectively, cosmic variance and shape noise. In the case of 3D underdensities, for clarity, we only show svf voids, but wvf and zobov voids show the same trends. We find that the S/N is boosted more when increasing the survey volume compared to when increasing the density of source galaxies by the same factor, although for troughs and svf_2d both strategies lead to similar S/N gains. In the case of 3D voids, shape noise represents a small fraction of the error budget, so having larger survey volumes brings the largest gain. For 2D underdensities, the shape noise error contribution is larger than in the case of their 3D counterparts, but it is still smaller than uncertainties due to sample variance and uncorrelated line-of-sight structures.

Figure 10.

Comparison of which sources of error dominate the S/N of the difference in tangential shear between F6 and GR. Each colour corresponds to a different void finder; for clarity we only show the svf results for 3D voids. The solid curves show the S/N for the lsst survey; the dashed lines show the change in S/N if instead the survey would have the same number of lensing sources but cover a four times larger volume; the dotted lines show the S/N if instead the survey would have four times more lensing source galaxies.

Open in new tabDownload slide

5 DISCUSSION AND CONCLUSIONS

We have carried out a detailed comparison of the extent to which voids can test models of modified gravity. Underdense regions are typically unscreened in modified gravity theories, which suggests that voids can be potentially good discriminators of such models. There are multiple ways of defining voids, with different algorithms designed to identify structures at different scales and with different geometries. Moreover, different void finders may have differing systematic and statistical error budgets. This motivated the need for a comprehensive comparison of different void-finding algorithms in order to precisely quantify which voids are best suited to test modified gravity models in light of the coming data from big galaxy surveys. In particular, we would like to know what level of constraints future surveys such as euclid and lsst can lead to when using void statistics. In this analysis, we compare six existing and new methods to identify cosmic voids, with three algorithms finding underdensities in the 3D galaxy distribution and three in the line-of-sight projected 2D galaxy distribution.

We expect that the constraining power of different void finders to be model dependent, so this work focuses on a class of very popular theoretical models: f(R) gravity proposed by Hu & Sawicki (2007), which is a representative example of the general models known as chameleons (Khoury & Weltman 2004). We studied two f(R) models, F5 and F6, with F5 corresponding to stronger deviations from GR than F6. The underlying philosophy is to use this very example to quantitatively understand the future constraints given by void statistics, and the conclusions from this study will be indicative for the more general models.

We improved upon previous studies, which used either dark matter particles or dark matter haloes as tracers, by identifying the voids using tracer galaxies, and therefore special care has been taken to create realistic mock galaxy catalogues. We have used a five-parameter HOD model to populate dark matter haloes with galaxies. The GR HOD uses the Manera et al. (2013) parameters, which reproduce the sdss cmass-sample galaxy clustering, while the f(R) HOD parameters were tuned to result in the same galaxy number density and projected two-point correlation functions as the GR ones. The resulting number densities of galaxies are n_g = 3 and 5 × 10⁻⁴ (Mpc/h)⁻³ for z = 0.5 and z = 0, respectively, and we leave a study of the effect of varying n_g (which requires higher-resolution simulations than used here) to future work.

The main conclusions of this work are as follows:

• 1. Void abundance The abundance of voids is sensitive to the tracers used to identify them. For example, Cai et al. (2015) found that the different f(R) models predict very different abundances for voids found using the dark matter field directly. For sparse tracers, such as dark matter haloes and galaxies whose number densities are in the region of ∼5 × 10⁻⁴(h⁻¹ Mpc)⁻³, we find the same void abundance across all models when matching the number densities, and the agreement becomes even better when we further match the galaxy correlation functions (see Fig. 4). Thus, void abundances do not have additional discrimination power of modified gravity models once the number density and the correlation function are the same across all models.

  The above conclusion holds when comparing voids identified using the same method for all models, and there can still be a lot of variation if we compare the abundances from different void finders.

• 2. Void profiles Similar to void abundances, the galaxy number density profiles of voids are almost identical across all models when we have matched their two-point correlation functions (see Fig. 5), and it suggests that the observed void galaxy number density profiles cannot be used to distinguish the models studied here.

  However, the modification of gravity does affect the distribution of the underlying dark matter field. The agreement of the galaxy correlation functions in the different models has been achieved by tuning the HOD parameters – which are empirical parameters describing how galaxies populate dark matter haloes from a simulation – in the models, and this means that these models must have different galaxy bias in order to match the same observational data. When studying the void matter density profiles, we find that voids in f(R) gravity are more underdense (see Fig. 6) due to the presence of a fifth force that evacuates underdense regions more efficiently.

  Gravitational lensing is a way to directly probe the total matter distribution between the source and the observer, and so we compared the stacked lensing signal (tangential shear) by voids between f(R) and GR gravity to find that the former generally predicts stronger void lensing (see Fig. 8). The tangential shear profiles of voids depend on the way in which the voids are identified. In particular, we confirm that among the 3D void finders, the svf gives rise to the strongest model difference in the tangential shear profile, since, by construction, spherical voids show a stronger density variation with distance at the void edge. We also find that 2D voids, in particular tunnels and 2D svf, show both stronger lensing signals and stronger model differences compared to 3D voids.

• 3. Constraining power of void lensing in future surveys In order to quantify the significance of the lensing signals for future galaxy surveys such as euclid and lsst, and to assess more accurately the potential of using voids to test chameleon-type models, we have performed a comprehensive calculation of the S/N of void tangential shear to distinguish f(R) and GR models. Our analysis has taken into account the major sources of uncertainties: the number densities of voids in the survey volume, the number density of source galaxies, line-of-sight projection effects, and the shape noise due to the intrinsic ellipticity of source galaxies.

  We find that 2D voids are the most promising underdensities for probing f(R) gravity. Of these, the tangential shear by tunnels has the largest constraining power, with a euclid-like survey being able to distinguish the f(R) gravity F5 and F6 models at a 50 and 7σ confidence level, respectively (see Fig. 9). The lsst data will have an even higher constraining power corresponding to a confidence level of 80 and 11σ for the F5 and F6 models. The 2D svf voids have somewhat less distinguishing power than tunnels, with troughs being the least constraining of the 2D voids.

  We also find that the weak lensing signal of 3D voids has a poor power to distinguishing f(R) models. The most promising of them, the svf, applied to the lsst survey can constrain the F5 and F6 models to a confidence level of 2 and 0.3σ. This distinguishing power is significantly lower than found in previous literature. The discrepancy with previous studies is due to (1) the use of the z = 0.5 snapshots, which show smaller differences between f(R) and GR than present at z = 0, and (2) our inclusion of additional error sources, such as shape noise and the contribution of uncorrelated line-of-sight large-scale structure.

  We investigated the dominant void lensing uncertainty sources for the lsst survey to find that while source galaxy shape noise is less important than errors due to sample variance and line-of-sight large-scale structures, it still has a significant contribution (see Fig. 10). Thus, the prospects of using void lensing to distinguish f(R) models from GR can be best boosted by increasing the survey volume. The best way of increasing the survey volume is by extending the sky coverage of the lsst survey; going to higher redshift will cover more volume, but may not help as much because at higher redshift the difference between models also decreases.

  This paper uses a conservative sample of tracer galaxies with a number density corresponding to that of the cmass galaxy sample. euclid and lsst are expected to have a larger number density of galaxies that will allow for the identification of more voids, but with smaller sizes. It remains to be studied how the larger density of tracers can affect the constraining power of void weak lensing. We leave this for future work, since it needs higher-volume simulations than used here with much better mass resolution, the latter needed to resolve the lower mass dark matter haloes that host galaxies fainter than in cmass. Another important question, which is not addressed in this paper, is the potential degeneracy between the effects of modified gravity and other cosmological parameters such as Ω_m and σ₈ (Cai et al. 2015). In f(R) gravity, the convergence power spectra show a scale-dependent enhancement compared to GR predictions (e.g. Tessore et al. 2015; Li & Shirasaki 2018), which suggests that the lensing signature of f(R) models is not degenerate with cosmological parameters, as explicitly checked by Shirasaki et al. (2017). Extending this conclusion to void lensing is non-trivial because the convergence power spectrum represents an average over the entire volume while voids sample mostly underdense regions, and, furthermore, there is an environmental dependence of the fifth force in f(R) models. Understanding the degeneracy between cosmological parameters and void lensing will involve new simulations to be carried out by varying these cosmological parameters and will be left for future work.

Acknowledgements

The authors would like to thank the anonymous referee for detailed comments that have helped us improve the paper and to Mark Neyrinck for insightful discussions. MC was supported by Science and Technology Facilities Council (STFC) [ST/L00075X/1, ST/P000451/1] and by an ERC Advanced Investigator grant COSMIWAY [GA 267291]. E.P. is supported by CONICYT-PCHA/Doctorado Nacional (2017-21170093). S.B. is indebted to the generous grant of an ITC Fellowship from Harvard University and, formerly, the STFC through grant ST/K501979/1. Y.C. was supported by the European Research Council under grant numbers 670193. N.P. received support from Fondecyt Regular 1150300, BASAL CATA PFB-06, Anillo ACT-1417. B.L. is supported by an European Research Council Starting Grant (ERC-StG-716532-PUNCA) and STFC (ST/L00075X/1, ST/P000541/1). This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement no. 734374. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1, and Durham University. DiRAC is part of the National E-Infrastructure. Part of the analysis was done on the Geryon cluster at the Centre for Astro-Engineering UC. The Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130008 provided funding for several improvements to the Geryon cluster.

Footnotes

REFERENCES

APPENDIX A: SELECTION OF TUNNELS AND TROUGHS

We are interested in studying the signature of f(R) modified gravity in underdense regions, but not all tunnels and troughs correspond to underdense regions. Here, we present one simple method, which can be applied to observations, of segregating tunnels and troughs according to their projected density contrast, Σ/〈Σ〉. For tunnels, their radius is correlated with their projected density, as shown in the top panel of Fig. A1. Large tunnels correspond to regions sparsely populated with galaxies and thus have low Σ values; in contrast small tunnels are found in regions crowded with galaxies, typically corresponding to high-mass haloes. A simple cut in radius, R_p eff ≥ 1 h⁻¹ Mpc, represents a good compromise between including as many objects as possible while still considering mainly underdense tunnels.

Figure A1.

The projected matter density, Σ, profile of tunnels (top) and troughs (bottom) at z = 0. The tunnels are split according to their effective projected radius, R_p eff. We use only tunnels with R_p eff ≥ 1 h⁻¹ Mpc that correspond to underdense regions. The troughs are split according to the number of tracer galaxies they contain. We use only troughs with N_gal ≤ 2 that correspond to a similar projected underdensity as the tunnels.

Open in new tabDownload slide

All troughs are defined to have the same radius, 2 h⁻¹ Mpc, so to discriminate between under- and overdense ones we use the number of galaxies inside each trough. Troughs that contain few galaxies are typically found in regions of low projected matter densities, as shown in the bottom panel of Fig. A1. The four bins in trough galaxy count were selected to contain roughly 20, 30, 30, and 20 per cent of the population, respectively. Our analysis uses only troughs that contain at most two galaxies inside them, which corresponds to 22 and 30 per cent of objects at z = 0 and z = 0.5, respectively. The z = 0.5 HOD has a lower galaxy number density and thus the same cut in trough galaxy count corresponds to a larger fraction of the population. The threshold of two galaxy counts in troughs is motivated by finding a compromise between large differences in the Σ profiles between f(R) and GR, which are maximal for the most underdense troughs, and including a sufficient fraction of the population.

APPENDIX B: LENSING COVARIANCE MATRIX

Figure B1.

The sample variance correlation matrices, R_ij, of the differential surface mass density, ΔΣ, for the six void finders used in this work.

Open in new tabDownload slide

The shape noise covariance was calculated by generating a spatially random catalogue of source galaxies whose number density and intrinsic ellipticity distribution matches that of the target lensing survey. Furthermore, the source galaxies are oriented randomly, which corresponds to a null tangential shear signal. We calculate the stacked mean galaxy ellipticity for each void catalogue by correlating the void distribution with the source galaxies. We obtain 500 bootstrap realizations that are then used to compute the shape noise covariance, which is shown in Fig. B2. For all void finders, we find that the diagonal terms are dominant. This is especially the case for tunnels and svf_2d voids, which show the smallest degree of overlap (see Fig. 3). The remaining void catalogues have various degrees of overlap that results in many off-diagonal elements being non-zero.

Figure B2.

The shape noise correlation matrices for three of the void finders used in this work. The correlation matrices shown here are top panel for svf voids (we checked that all 3D voids have a similar matrix), central panel for tunnel (we checked that the svf_2d result looks almost identical), and bottom panel for troughs.

Open in new tabDownload slide