© ESO 2019

1. Introduction

Based on Einstein’s General Relativity (GR) theory, most of the observational features of the Universe on cosmological scales are nicely reproduced by the so-called cosmic concordance model or Λ cold dark matter (CDM) (Perlmutter 2003; Eisenstein et al. 2005; Kowalski et al. 2008; Rapetti et al. 2008; Hinshaw et al. 2013; Planck Collaboration VI 2019). Nevertheless, one of the main problems in modern cosmology is understanding the nature of dark energy. This is the exotic energy component with negative pressure that is required by the standard ΛCDM model to reproduce the accelerated expansion of the Universe, which is observed at low redshift (Riess et al. 1998; Perlmutter et al. 1998). The simplest candidate for dark energy is the cosmological constant, denoted by Λ, which is a free geometrical parameter of Einstein’s theory.

Some alternatives to the cosmological constant are the quintessence or dark sector interaction models, for instance, and running the vacuum expectation value. For a review, see Yoo & Watanabe (2012), Joyce et al. (2015), and references therein. We are here interested in another set of theories, that is, theories of modified gravity (MG). These aim to model the accelerated expansion of the late Universe by going beyond GR. Viable MG theories must display a screening mechanism, which consists of the weakening of the gravity modifications within dense environments, such as our Solar System, where GR has been exhaustively tested. Different families of MG theories are characterized by different ways to accomplish the screening effect (see, e.g., Brax et al. 2012; Koyama 2016).

When efficient screening mechanisms are available, MG theories are virtually indistinguishable from GR inside massive halos. On the other hand, long-range forces as well as cumulative effects due to different time-evolution paths can modify the large-scale spatial distribution and the abundance of halos. These observables have recently been studied in different contexts (e.g., Schmidt et al. 2009; Winther et al. 2015; Koyama 2016). In contrast, we here focus on a complementary scenario where screening mechanisms could be weak, leading to a detectable signature of modified gravity. Our goal is to test MG through the analysis of cosmic voids, which are the intermediate-sized and underdense regions that are left behind by the hierarchical structure formation of dark matter halos.

By definition, voids are underdense regions naturally bounded by an overdense or a background-dense wall. Methods to detect them include free-parameter geometrical algorithms intended to find non-spherical voids (Lavaux & Wandelt 2012; Neyrinck 2008; Jones et al. 2007). Spherical void finders are based on theoretical density thresholds (Colberg et al. 2008).

Many cosmological probes based on void analysis have been proposed in the past decade in the context of GR, including a number of observables such as ellipticity (Park & Lee 2007; Shoji & Lee 2012), density profile (Nadathur et al. 2014, 2019; Sánchez et al. 2017), and gravitational lensing (Sánchez et al. 2017; Krause et al. 2012). Moreover, void properties depend strongly on the void-finding algorithm (Jones et al. 2007; Chan et al. 2014; Elyiv et al. 2015; Frenk et al. 2016; Hamaus et al. 2017). Despite this dependence on the void-finding algorithm, the void analysis can be performed in a self-consistent way by calibrating the methods on mock catalogs, which plays the role of the theory, and this has been successfully applied to photometric samples (Nadathur et al. 2019) that improved the Baryonic Acoustic Oscillations (BAO) constraints. This type of analysis shows great promise if it were applied to future experiments, including DESI (Aghamousa et al. 2016) and Euclid (Scaramella et al. 2014).

The study of voids in the context of modified gravity has recently gained strength. This includes mainly the effect of MG on the void abundance (Li et al. 2012; Jennings et al. 2013; Lam et al. 2015; Cai et al. 2015; Voivodic et al. 2017) and in the lensing signature around voids (Barreira et al. 2015; Winther & Ferreira 2015; Baker et al. 2018). In these scenarios, the galaxy distribution is modified by the fifth force associated with MG, and in turn, this also affects the shape and abundance of voids (Zivick et al. 2015). MG theories must reproduce the Newtonian gravitational force within dense regions while displaying the characteristic repulsive effect of the cosmological constant in background density environments. The transition between these two asymptotic behaviors is a potential scenario to probe MG models. We therefore focus here on void-related observables, and compare the effects of different screening mechanism effects on the ΛCDM outcome.

Specifically, we aim to distinguish between the ΛCDM model and two MG theories. The first theory is the Hu–Sawicki f(R) theory (Hu & Sawicki 2007), which implements the chameleon screening mechanism. The second theory is the symmetron model (Olive & Pospelov 2008; Hinterbichler & Khoury 2010; Hinterbichler et al. 2011), which implements the screening mechanism that bears its name.

This study was carried out in the context of an idealized scenario, at z = 0 and with a fixed set of cosmological parameters whose variation would lead to already observed degeneracies. Generalizations, like those cited above, will be part of our future efforts.

The paper is organized as follows. In Sect. 2 we give an overview of the MG models we are interested in. In Sect. 3 we describe the N-body simulations we analyze in this work, as well as the void-finder algorithm we employ. In Sect. 4 we describe the three void-related observables we choose to analyze, how we estimate them from the N-body simulations, and the phenomenological fitting expressions we use to describe them. In Sects. 5 and 6 we show how well we can distinguish among GR and the two MG theories based on the abundance, linear bias, and density profile of voids. Finally, in Sect. 7 we present our conclusions.

2. Gravity models

2.1. Symmetron

The symmetron cosmological model (Hinterbichler & Khoury 2010; Hinterbichler et al. 2011; Olive & Pospelov 2008) is a scalar–tensor theory for a single scalar field ϕ. The action for a scalar–tensor theory can be written as

(1)

where S_m is the action associated with the standard matter fields ψ_i. The fields ψ_i are coupled to the scalar field ϕ through the Jordan frame metric , which is related to the Einstein frame metric g_μν by . For the particular case of the symmetron theory, we have

(2)

and

(3)

In this case, both A(ϕ) and V(ϕ) are symmetric upon the transformation ϕ → −ϕ. Here, μ and M are mass scales, λ is a dimensionless coupling constant, is the Planck mass scale, and g is the determinant of g_μν.

The variation of the action in Eq. (1), with respect to the Jordan frame metric provides the dynamical equation for the scalar field ϕ:

(4)

The effective potential V_eff, Eq. (4), has a minimum at ϕ = 0 in high-density environments, that is, where the local matter density ρ_m ≫ μ²M² ≡ ρ_SSB. In this case, V_eff holds the ϕ → −ϕ symmetry. On the other hand, in low-density environments, ρ_m ≪ ρ_SSB, the effective potential displays two minima at , breaking the symmetry of the model. Here is the expected value of ϕ for ρ_m = 0.

The free parameters of the symmetron model {μ,  M,  λ} are usually exchanged by the physical parameters associated with the scalar field for ρ_m = 0 (Winther et al. 2012):

(5)

which correspond to the range of the scalar field in Mpc h⁻¹ units (λ₀), the dimensionless coupling strength to matter (β), and the redshift for which the symmetry breaking occurs at the background level (z_SSB), respectively. z_SSB is related to the symmetron critical density ρ_SSB, for which the symmetry breaking takes place through the expression . Here is the matter background density at redshift zero.

2.2. f (R) gravity

The Hu–Sawicki f(R) model (Hu & Sawicki 2007) was first formulated in the Jordan frame in terms of the action

(6)

where ℒ_m is the Lagrangian describing the matter fields, and

(7)

In this case, , where H₀ and Ω_m0 are the Hubble constant and matter density relative to the critical value today. In order to recover the effective cosmological constant in the large curvature regime, it must be set c₁/c₂ = 6Ω_Λ/Ω_m.

This f(R) theory can be transformed into a scalar–tensor theory upon both the identification

(8)

and the frame transformation

(9)

In this case, the Compton wavelength or range of propagation of the scalar field at redshift zero is given by

(10)

where f_R0 can be expressed as a function of {c₂,  n} as

(11)

3. Methods

3.1. N-body simulations

The cosmological N-body simulations analyzed in this work were run with the ISIS code (Llinares et al. 2014; Llinares & Mota 2014), which is a modification of the GR RAMSES code (Teyssier 2002) to include MG models, with 512³ dark matter particles in a (256 Mpc h⁻¹)³ cubic box. The initial conditions correspond to a flat ΛCDM cosmology with parameters {Ω_c,  Ω_b,  h,  σ₈,  n_s} = {0.222,  0.045,  0.719,  0.8,  1} and without neutrinos. Massive neutrinos and modified gravity degeneracy has been pointed out and analyzed before (see Hagstotz et al. 2019; Schuster et al. 2019; Kreisch et al. 2019), but we here prefer to exclude massive neutrinos to single out the modified gravity effects.

We intend to analyze void properties within the framework of the symmetron and f(R) gravity. N-body simulations are highly computationally expensive, therefore it is not feasible to cover the full parameter space of these theories. On the other hand, we aim to show the MG effects on a recently explored void probe as clearly as possible, therefore we do not restrict our analysis to tighter constraints (Burrage et al. 2019; Tsujikawa 2008; Pogosian & Silvestri 2010). The MG cases we cover here are described by the parameters n = 1 and |f_R0| = {10⁻⁴,  10⁻⁵,  10⁻⁶} for the f(R) theory, and by β = 1, λ₀ = 1 Mpc h⁻¹ and z_SSB = {1,  2,  3} for the symmetron theory (see Tables 1 and 2). The set of simulations is complemented by the ΛCDM or GR case, where no gravity modifications are included. We assume here that the ΛCDM case corresponds to both the |f_R0| = 0 and the z_SSB = 0 scenarios.

Even though the effect of f(R) and symmetron gravities may change over redshift, the effect is expected to be stronger at z = 0. In particular, symmetron and GR must be indistinguishable for z >  z_SSB, while in linear theory the enhancement of gravity due to f(R) weakens for higher redshift (Schmidt et al. 2009). We here only analyze the z = 0 outputs and leave the redshift analysis for future works.

3.2. Void-finding algorithm

Non-spherical voids, such as those found by methods based on Voronoi tessellation, display a dependency of the inner density on void size. In these cases, the smaller the void, the emptier (Hamaus et al. 2014a), in contrast to the spherical model prediction, where each void has the same mean density of about 0.2 times the background density. Because denser voids can naturally be larger, they leave less room for small and emptier voids. The abundance functions of voids that are detected through Voronoi tessellation and fixed density are therefore considerably different.

On the other hand, the ellipticity distribution of non-spherical voids also depends on the void size (Park & Lee 2007; Shoji & Lee 2012). Using the non-spherical version of our void-finder algorithm, which we describe in the next paragraph, we found very little information about MG effects in the ellipticity distributions. Because of this, and because we are interested in the spherically averaged density profile of voids, we considered spherical voids to be more suitable for this analysis. In principle, this also helps us to parameterize the void density profiles more simply.

According to spherical expansion theory, in an Einstein-de Siter (EdS) cosmology, voids are predicted to have an average overdensity Δ_v = 0.2 that reaches as far as the void radius (see Appendix A of Jennings et al. 2013). Even though Δ_v depends on cosmology and MG parameters, we set it to its EdS value as is commonly done in the case of the density contrast for defining halos in simulations. Moreover, in a real application, we would not know the correct gravitational theory and cosmology to compute Δ_v a priori.

The void-finding algorithm starts by computing the dark matter density on a regular grid of size l. This is done by applying the cloud-in-cell (CIC) algorithm to the dark matter particles. The overall process consists of three stages that we describe below.

(1) Initial grid spheres. First, a sphere centered on each grid cell j is grown until radius r_j, where the mean density of the sphere reaches the critical value , which is the expected density for spherical voids at redshift z = 0 (see Jennings et al. 2013). Here, is the mean density of the dark matter particles in the cosmological box. We refer to these spheres as grid spheres, and some of them are refined in the next step.

(2) Adaptive refinement. The second stage improves the initial estimate of the radius and the center position of the grid spheres. This improvement is accomplished using the particle positions instead of the density of the grid cells as was done in the first stage to compute the densities. This step aims to maximize the size of the sphere that is refined and makes the result robust regardless of the first guess given by the first stage.

We improved upon stage one by (i) growing spheres with an average density with a critical value, not only at the current position x_j (initially the center of the cell), but also at the corners of a cube of side d around x_j. (ii) When one of the corners maximizes the size of the sphere that is refined, we moved x_j to that position. Otherwise we did not update x_j, but instead reduced d to half of its current value (the initial value for d is the grid side l). For a given sphere j, we iterated this refinement process (steps i and ii) until d reached the threshold defined as the minimum value between 0.125 Mpc h⁻¹ and 1% of the current sphere radius r_j.

When the radius r_j and position x_j of a given sphere j were refined, we set r_k = 0 for every grid sphere k whose center x_k was closer than 0.9r_k to x_j. This was done to avoid duplicates because under refinement, these grid spheres k would converge to the same values r_j and x_j. We stopped refining grid spheres when none of them had a radius larger than 2l. After this stage, we had what we call the catalog of void candidates (grid spheres whose radius and position were refined).

(3) Family casting. We gathered void candidates (denoted here by C) into families (which we call voids and denote by V) by applying the linking procedure described below. In this way, a candidate added to a family becomes a void member, while a void is identified as the collection of its members.

First, we assigned the largest void-candidate to the first family. Then we considered the next largest candidate C, whose core was defined as the sphere around its center with a radius that is 70% of the candidate’s radius. Casting then proceeded as follows:

– If none of the already identified voids overlapped the core of the candidate (meaning that C is isolated enough from any already detected void), we created a new family and assigned C to it.

– If the center of C was inside an already identified void V and no other void overlapped the core of C, we assigned the candidate C to the void V.

– Otherwise, we removed C from the candidate catalog because either it was unclear to which of the already detected voids the candidate C belongs, or because C was not isolated enough to define a new void. If C was not discarded, it could become the linking piece of a bilobe (dumbbell-shaped) void, which is not the type of void we are interested in.

After all the spheres in the candidate catalog were cast, roughly one-third of them were discarded. The remaining two-thirds of the original candidates were gathered into families (which we call non-spherical voids) with an average of about three spheres per void. The largest voids have nearly one hundred spherical members, while many small voids have only one spherical member (which is a consequence of the finite number of CDM particles in the simulation).

In Fig. 1 we illustrate the family-casting procedure. The condition for turning a candidate C into a new void is that none of the already detected voids touches the core of C. This is the case of c_r, for which a new family is created, turning c_r into the first member of this new family (void). The two conditions for assigning a candidate C to an already detected family V are that the center of C has to be inside V, and the overlap between the core of C and any other void has to be zero. This is the case for c_m, c_n, and c_t. In this case, c_m and c_n are assigned to family V_i during the casting process, while c_t is assigned to family V_k. The center of c_t is inside a single void (family), and no other family overlaps the blue core of c_t. Finally, candidates {c_p,  c_q,  c_s,  c_u } are discarded during the casting process. The centers of c_s and c_u are not inside any void, but they are not isolated enough from already detected voids (V_k touches their cores). On the other hand, c_q is not isolated enough for becoming a new void, and its core is touched by both V_j and V_l, therefore it is unclear to which family c_q would belong. Likewise, the core of c_p has reasonable overlap with both V_j and V_k.

Fig. 1. Illustration in 2D of the assignment of void candidates to families (voids). We show voids that have already been identified {Vi,  Vj,  Vk,  Vl}, with both Vk and Vl having more than one member. The centers of each member in a family are inside the radius of another member belonging to a previous generation of the same family (we highlight this by showing the center of the family members, except for the very first detected member of each family). Blue circles denote void candidates and blue shades represent their cores (which correspond to 70% of the candidate’s radius).

3.3. Non-spherical and spherical voids

The non-spherical void catalog consists of all the families that were founded by applying steps (1)–(3) described in the last subsection. As a direct consequence of the finite number of CDM particles in the simulation, the resulting non-spherical void catalog can be seen as being composed of a subsample of small voids that are dominated by spherical members and a subsample of large voids that are dominated by non-spherical members.

We extracted a spherical-void catalog from the non-spherical catalog by defining a spherical void as the largest member of each non-spherical void. In this case, each spherical void has a mean density equal to by construction, and its radius is denoted by r_0.2.

Extracting a spherical-void catalog from a non-spherical one instead of directly detecting spherical voids from the very beginning weakens the problem of breaking highly non-spherical voids into several spherical pieces, which causes miscounting. Furthermore, in this case, the additional voids would inevitably be bounded by other voids, which would create a subsample with different features than those proper of voids that are bounded by overdensities. In order to work with a purer sample, we chose to analyze the catalog of spherical voids alone.

4. Void observables

In this section we present the phenomenological expressions we used to describe the void abundance f(σ), the void density profile ρ(r,  r_0.2), and the void-matter linear bias b(σ). Here, σ² = σ²(R) is the variance of the linear matter power spectrum smoothed on a scale R, which was computed for the ΛCDM and MG cases as described in Voivodic et al. (2017).

The expressions in this section were chosen with the aim of describing the voids properties (just as the Navarro-Frenk-White halo density profile, Navarro et al. (1997), or some halo mass functions, see, e.g., Tinker et al. (2008), which are very useful fits for describing N-body simulations) with enough accuracy in order to constrain the MG parameters from simulated data, which is the main goal of this work. They can lack a direct origin from first principles, but they are inspired in some theoretical results like the peak background split for describing the void-matter bias, or the excursion set formalism for the void abundance.

The specific expressions for the different void properties already depend on cosmology and gravity through the σ function, but we found that the free parameters in these expressions depend on MG. We take this MG dependence into account through the functional form

(12)

which are slight modifications of linear functions in log₁₀|f_R0| or z_SSB. If they were pure linear functions (γ₃ = γ₆ = 0), the quality of the fit for the z_SSB = 0 case would be reduced, and our fiducial value for |f_R0| would be ∼10⁻⁸ instead of 0. On the other hand, γ₄ could be written as γ₄ = γ₁ + γ₂log₁₀(γ₃) in order to ensure a unique description for the GR limit for both f(R) and symmetron, but we do not need to do this because we analyzed both theories independently.

The functional for in Eq. (12) appears recurrently in this section. We therefore include the MG dependence in an almost linear way because all the MG cases we study here are quite close to GR. Mainly because no 1-to-1 map for f(R) and Symmetron is available, we were unable to find a common parameterization for the MG dependence. Furthermore, this non-equivalence between these theories is the main motivation for us to search for a probe to distinguish between them.

4.1. Void abundance

It has been shown (e.g., Voivodic et al. 2017) that the abundance of spherical voids is well described by the excursion set formalism and parameters δ_v and δ_c that come from spherical expansion or collapse theories. The voids of Voivodic et al. were grown centered on positions given by particle coordinates of the minimum of the density field given by the Voronoi volume of each particle.

We defined the center of the voids in order to maximize the void radius, therefore we obtained voids that are larger than those described in Voivodic et al. (2017). Moreover, because we did not set the center of our voids on the Voronoi volume particles, we detected more small voids than Voivodic et al. (2017). As a result, our void abundance is not well described by the excursion set model described in Voivodic et al. (2017). Instead, we describe the void abundance with a phenomenological formula, using the same functional form of the halo mass function as in Tinker et al. (2008),

(13)

where need to be fit. After the fitting process, we found that A, b, and c can be considered common constants for the theories we analyzed. On the other hand, we needed to introduce an explicit dependency on the gravitational theory through , Eq. (12). This allowed us to account for more large voids when MG is stronger (due to a higher level of clustering), therefore leaving less room to small voids.

On the other hand, our void-detection algorithm is not sensitive to the void-in-cloud process, which prevents small voids from surviving in dense regions (Jennings et al. 2013). From our point of view, the void-in-cloud process is more suitable for a continuous description of the dark matter density field and is not a property of a discrete system such as an N-body simulation or a galaxy catalog.

In a discrete system, our algorithm will always find more voids when we scan over smaller scales (as long as there are enough tracers available). This outcome is not compatible with the theoretical predictions that take the void-in-cloud process into account, but instead, it is more similar to the halo abundance functional form, Eq. (13). Furthermore, the simplest model for the void abundance corresponds to the Press-Schechter result, which has the same functional form for both voids and halos. Therefore we expect that a more elaborate halo abundance functional form is required to describe the void number-counts in the context of our analysis.

In Fig. A.5 we show the measured void abundance and the best fit of Eq. (13) for the ΛCDM, f(R), and symmetron N-body simulations. The errors on the measurements, estimated as the variance of the octant-subsamples in the simulated box, closely follow the Poisson noise expectation of the entire sample. We therefore used Poisson errors for the void abundance.

4.2. Void density profile

The void density profile was estimated as the mean of stacked voids traced by the dark matter particles as

(14)

with

(15)

where m and r_i are the mass and position of the dark matter particles, while 2δr = 0.05r_0.2 is the thickness of the shells we used to sample the density profile.

We split the void catalog by void sizes into seven intervals of r_0.2. Then, we rescaled the individual density profiles from ρ(r) → ρ(r/r_0.2). Finally, we stacked every rescaled density profile belonging to the same size interval. The errors on the measurements were estimated from the variance of the octants in the simulated box.

We used the following phenomenological expression (Hamaus et al. 2014a) to describe the void density profile:

(16)

where is the background dark matter density, δ₀ is a constant, and the parameters {α,  β,  c,  G} depend on the void size as well as on the free parameters of the MG theory. We parameterized them as functions of the MG parameter x and σ(r_0.2) as follows:

(17)

and have the same functional form as in Eq. (12) in the f(R) case, while they are linear functions of (x = z_SSB) in the case of the symmetron theory. All the subindexed coefficients are constants, but they have slightly different values for f(R) and symmetron because we analyzed each theory independently in this work. In Eq. (17), σ ≡ σ(r_0.2) contains cosmological and MG information as well as the void size, while x is the MG parameter itself (see Eq. (12)). As a result, the f(R) (symmetron) theory has a total of 11 (13) parameters to be fitted.

We fit for the free parameters in seven stacked profiles for each one of the different values of |f_R0| or z_SSB. The seven stacked profiles for each MG case include all the voids in the same bin. These seven bins share the same length in log(r_0.2), and they split the void sample r_0.2 ∈ (1.5 Mpc h⁻¹,  14 Mpc h⁻¹) by size. The value of the radius reported for a given stack corresponds to the mean of the void radius in each bin. On the other hand, the density was measured for {150, 140, 130, 120, 110, 100, 90, 80} spherical shells of thickness 0.11/r_0.2 for the smallest to the largest voids associated with the seven bins defined above. In the next paragraphs, we describe the heuristic choices we made for the particular parametrization in Eq. (17).

(α, β) For r ≳ 2 r_0.2 the Expression (16) falls like r^{α − β}, in our case as r^{(1 − β₀)α} with (1 − β₀) < 0 and α >  0, describing the tail of the profile. This tail is steeper for larger voids because it is the linear void-matter correlation function for larger scales, and we approximated this effect with the σ dependence of α. The MG dependence in the term was added because the steepness and amplitude of the matter-void correlation function is systematically higher for stronger modifications in the gravity force.

(c) The value of c^−β sets the position of the peak in Expression (16). Our voids are defined to have the same mean density (below the background level) at r = r_0.2, which causes their peaks to stand away from the center in the cases of higher walls (smaller voids or stronger MG). This effect is taken into account by letting c^−β increase with σ (decrease with r_0.2).

(G) Finally, the second term in G allows us to move from a positive to a negative density around and beyond the wall of the void when we consider larger voids. Again, lets us mimic the effect that MG shows stronger clustering on small scales, which causes the change of sign to occur for higher values of r_0.2.

Some void density profiles are shown in Figs. A.1 and A.2. We show for three different void sizes and the different MG cases.

We set the density at the void center δ₀ to a constant value because we did not take the inner part of the profile into account in this analysis. We computed δ₀ as being the average, over all the void sizes and gravity cases, of the central density by fitting a power law in r/r_0.2 plus a constant to the inner void profiles. This was done because Eq. (16) cannot describe the inner and outer parts of the profiles very well simultaneously for the spherical voids because there is an abrupt change in the measured profiles around r = r_0.2 that is due mainly to the steepness of the walls. Still, this functional form follows the shape of density profile (of the voids analyzed in this work) better than other proposals in the literature (van de Weygaert & van Kampen 1993; Maggiore & Riotto 2010; Colberg et al. 2005; Lavaux & Wandelt 2012; Bolejko et al. 2013; Ricciardelli et al. 2013, 2014).

4.3. Matter-void linear bias

The linear bias between the dark matter density field and the void distribution on large scales was estimated as

(18)

where P_mm is the matter power spectrum and P_mv is the matter-void cross spectrum. An alternative definition for the void bias would be . We did not use the last expression for two reasons. First, the estimation of P_vv is prone to shot-noise much higher than P_vm. Second, on large scales, P_vm changes sign for voids with r_0.2 ∼ 7 Mpc h⁻¹, while P_vv is always positive by definition. This change of sign in b is expected by the peak-background-split (PBS) prediction (Chan et al. 2014) and has been seen in simulations (Hamaus et al. 2014b). Therefore, using P_vm to estimate the void-matter bias allows us to use a wider range of void sizes, which gives better constraints.

The PBS approach (though mainly for describing halos) consists of splitting the total matter density field into two independent components: a long-wavelength background contribution, and a short-wavelength peak contribution. In this scenario, the number of peaks will change differentially with the background contribution for a given density threshold that defines the peaks. This difference can be translated into a linear bias between the total density field and the number density of the tracer (the peaks).

When a universal halo mass function there is available, the PBS can be used to compute a self-consistent linear halo-matter bias. A mass function as predicted by the excursion set formalism and a PBS bias both describe roughly simulated data and provide indications for constructing more accurate phenomenological counterparts. The same methods have been broadly applied in the literature to the void description (see, e.g., Sheth & Van De Weygaert 2004; Chan et al. 2014), which we take as a frame to set up more flexible expressions to tightly fit the N-body simulations including the MG effects.

We parameterize the linear void-matter bias as

(19)

where a and d are constants, while has the same functional form as in Eq. (12) for the f(R) case, but is a linear function of (x = z_SSB) for the symmetron case.

The results of the next section show that MG constraints based on the linear bias are weaker than those from the abundance and the density profile analyses because the length of the simulation box is 256 Mpc h⁻¹, giving a minimum Fourier mode k = 0.025 h/Mpc and then a poor sampling of large scales. We used the first five linear bins of length 0.05 Mpc h⁻¹ to fit the linear trend of the matter-void bias, as shown in Fig. A.4 for the ΛCDM case. The errors associated with spectra P_vm and P_mm were estimated from the variance of all the modes that contribute to a given bin in Fourier space, while the errors on the linear bias come from the fit of Eq. (18) as a linear function of k² for k <  0.25 Mpc h⁻¹ (see Fig. A.4).

5. Constraining modified gravity

The first part of our analysis consisted of using the simulations to fit for the free coefficients of the phenomenological models for the void abundance Eq. (13), density profile Eq. (16), and void-matter linear bias Eq. (19), assuming the fiducial values for the MG free parameter in each case, see Tables 1 and 2. Then, we applied the derived phenomenological models to the same simulations in order to see how well we can recover the free MG parameter in each case. Figures 4 and 5 and Tables 3–10 summarize the results.

We highlight that the constraints associated with the linear bias analysis are considerably weaker than those associated with the abundance or with the density profile analyses. For this reason, we present these result in separate plots, see Figs. 2 and 3, for example. The bias analysis shows that we can distinguish between GR and MG with parameters larger than |f_R0|≳10⁻⁶ or z_SSB ≳ 1 for the f(R) or symmetron scenarios, respectively. The results of the abundance, linear bias, and density profile analyses are compatible with each other and also with the fiducial values within two standard deviations.

Fig. 2. Constraints on the f(R) parameter fR0 derived by the matter-void linear bias analysis.

Fig. 3. Constraints on the symmetron parameter zSSB derived by the matter-void linear bias analysis.

Fig. 4. Posterior distributions for the free parameter of the f(R) theory when we analyzed the {ΛCDM,  f6,  f5,  f4} simulations from top to bottom. We note that |fR0| constraints from abundance and density profile are compatible at {1,  2,  1,  1} σ for the {ΛCDM,  f6,  f5,  f4} cases; see Table 12.

Fig. 5. Posterior distributions for the free parameter of the symmetron theory. The zSSB constraints from abundance and density profile are compatible at {1,  1,  2,  2} σ for the {ΛCDM,  A,  B,  D} cases, respectively, see Table 12.

In Appendix A we show the best fits of the joint analysis for the different simulations. In these plots we show that the symmetron and the f(R) effects over the three observables point in the same direction:

– Stronger MG means a stronger gravitational force, which results in more clustering in intermediate scales (gravity inside halos would be the same). Accordingly, it also means more high-mass halos because the more clustered they are, the more merging occurs. This boost in halo clustering would translate into large voids becoming larger in MG. As a consequence, there is less room for small voids, as we show in Fig. A.5.

– As was shown by Voivodic et al. (2017), the linear power spectrum is larger in MG at small scales. This effect translates into a matter correlation function that is also larger at intermediate and small scales in stronger MG cases. Therefore the walls of the voids will be higher in MG gravity, as we show in the simulations, see Fig. A.2. In consequence, because all our voids have the same mean density, the wall of the voids would also be steeper for stronger MG cases.

– When the linear void-matter bias is considered to depend on the void abundance, the natural consequence of having more large voids (and fewer small voids) in MG than in GR would be a higher bias. The steeper the abundance function of a tracer, the more positive the tracer-matter bias in the PBS approximation.

These features suggest that the symmetron and f(R) theories could be indistinguishable from each other when we consider only the bias, density profile, and abundance of void analyses. To address this question, we applied the symmetron analysis to the f(R) simulations and vice versa. The results are shown in Sect. 6.

6. Distinguishing among gravity theories

In order to assess the level to which the three gravity scenarios GR, f(R), and symmetron can be distinguished, we applied the f(R) analysis to the symmetron simulations and vice versa. This can show how well MG constraints derived from a correct model compare to those from an incorrect theory choice.

The results are summarized in Figs. 6 and 7. The constraints from the abundance, density profile, and bias probes are less consistent with each other when we perform the analysis using the incorrect MG model (see Table 12). Similarly, as shown in Table 11, the values of χ²/d.o.f. are higher when we assume an incorrect model.

Fig. 6. Recovered value for |fR0| for the symmetron simulations.

Fig. 7. Recovered value for zSSB for the f(R) simulations.

6.1. Analyzing symmetron simulations using f(R) theory

In the case of a joint analysis using all void probes, the lowest χ²/d.o.f. are {1.10,  1.78,  2.53} times higher for the {A,  B,  D} simulations when we apply the f(R) analysis than when the correct theory is used (see Table 11). Additionally, the best-fit values for |f_R0| from the abundance and density profile analyses disagree with each other at more than {5,  8,  4} standard deviations for the {A,  B,  D} cases (see Table 12). The same result is shown in Fig. 6.

6.2. Analyzing f(R) simulations using symmetron theory

In this case, the lowest χ²/d.o.f. from the joint analysis is {1.08,  1.31,  1.89} times higher for the {f6,  f5,  f4} simulations when they are analyzed with the symmetron theory than when f(R) is used (see Table 11). Likewise, the mean best-fit values for z_SSB from the abundance and density profile analysis disagree with each other at more than {6,  9,  13} standard deviations for the {f6,  f5,  f4} cases, as shown in Table 12 and Fig. 7.

6.3. Discussion

We recall that when the simulations are analyzed with the correct theory, the abundance and density profile analyses agree within two standard deviations, as discussed in Sect. 5. The results involving χ²/d.o.f. from the previous subsections show that we cannot reasonably distinguish between symmetron and f(R) modified theories for the weaker MG cases we considered (i.e., z_SSB = 1 and |f_R0| = 10⁻⁶). The values of χ²/d.o.f. are only marginally higher when an incorrect MG model is used.

One interesting point is that applying the correct model provides more consistency between the density profile and the abundance tests, however. This may be due to the fact that different screening mechanisms affect the void abundance and the profile differently. Therefore, the modeling from an incorrect MG model may effectively describe the simulated data for each void probe, but leads to conflicting values for the best fits. In a real data analysis, this difference might indicate that an incorrect gravity model is used.

Differences might also point to observational systematics that affect multiple observables differently. In a real analysis, we detect voids from discrete galaxies, not the dark matter field. Therefore we expect real void catalogs to suffer from a number of observational effects that make their use for cosmological purposes much harder than what we infer here. These effects can be collectively cast in the so-called void selection function, described by completeness and purity functions, for example, which are specific to each void-finder algorithm and survey data. See Aguena & Lima (2018) for a discussion of how imperfect knowledge of the cluster selection function, for instance, affects cosmological constraints derived from galaxy clusters. We expect similar effects for voids.

The key question is how well we can simultaneously constrain the parameters of the void selection function along with cosmological parameters of interest. A particular void finder may fail to detect true voids of a given void size (lowering completeness) and may also produce false voids (lowering purity). Misidentifications are more likely to occur for small voids or voids whose walls contain halos with a small number of galaxies. For sufficiently low-mass halos, the average number of galaxies per halo becomes ∼1 or lower. In this limit, regions devoid of galaxies are not necessarily empty of dark matter. As a result, there is no 1:1 correspondence between dark matter voids and galaxy voids. In order to account for these effects in simulations, we need to populate dark matter simulations with galaxies and replicate all observational features such as survey mask, flux, or magnitude limit cuts. This is beyond the scope of this work, but is necessary for accurate cosmological constraints from realistic void surveys.

For the cases where MG effects are stronger (|f_R0| = 10⁻⁴ and z_SSB = 3), applying the correct MG theory analysis provides much better consistency between the different tests, and the minimum χ²/d.o.f. is considerably lower (compared to the incorrect MG theory). Therefore at the level of N-boby simulations, we can safely distinguish between symmetron cases with z_SSB ≳ 2 and f(R) cases with |f_R0|≳10⁻⁵.

Interestingly, Figs. 3–7 show that any of the three void probes can tell us whether gravity is modified or not because none of the MG cases is compatible with ΛCDM, even when an incorrect MG theory is assumed. Moreover, when we analyze the ΛCDM case, we recover |f_R0| = 0 or z_SSB = 0. This means we will be able to know whether the Universe is ruled by GR or MG after applying any of the void analyses we considered here, even though it may be challenging to distinguish between the weakest cases of f(R) and symmetron by the joint analysis of the linear bias, density profile, and abundance of voids alone. On the other hand, if MG is stronger than the weakest cases considered here, the individual analyses may indicate differences in their best fits, forcing us to consider other MG models.

Finally, we note that the degeneracy between the weakest cases of MG analyzed here might in principle be broken by combining our void analysis with information from halo properties, such as their abundance, bias, and profiles. Recently, the splashback radius (Adhikari et al. 2018; Contigiani et al. 2019) and the turnaround radius (Lopes et al. 2019, 2018) of halos have been shown to be signficantly affected by MG effects.

7. Conclusions

We used a spherical-void finding algorithm to construct void catalogs in N-body simulations of GR as well as f(R) and symmetron theories. We measured the abundance, bias, and profiles of these voids and modeled our measurements with phenomenological fitting formulae. We then used these expressions and the simulated data to assess how well MG models can be constrained from void properties.

The void-finding algorithm we used can be described by the following main points: First, we searched for underdense spheres with a mean density given by (which defines the void radius r_0.2), centered on the cubic pixels of a regular grid. The pixel size was set here to be half of the mean particle distance, therefore we can take advantage of the full resolution of the simulations. Next, we maximized the radius of each sphere by refining the position of its center. Non-spherical voids were defined as the union of spheres with a neighbor-to-neighbor overlap larger than a given threshold, while spherical voids were selected as the largest sphere of the non-spherical ones.

The voids we found show an abrupt change in the density profile around the void radius. This makes it difficult to fit both the inner and outer regions of the void profile simultaneously. On the other hand, the height of the void walls increases in MG scenarios, and its variation with MG parameters is much more significant than the inner profile variation. Therefore, we chose to use only the outer part of the profile in our analysis.

Clearly, the changes due to MG observed in void properties are connected to those observed in halo properties. It is well known that the matter power spectrum and the halo properties change quite significantly as a function of MG parameters (e.g., Schmidt et al. 2009; Wyman et al. 2013). Viable MG models increase gravity effects, causing massive halos to become more abundant and voids to become emptier. Because the most underdense regions (δ ∼ −1) in a GR scenario cannot be much emptier in MG, the inner regions of voids do not change significantly, and most of the void profile modifications are concentrated around the void walls. Likewise, void radii are larger in MG than in GR, which is also compatible with the MG effects on halo properties. The more massive halos are more clustered on void walls in MG, producing higher walls and larger voids.

We parameterized the void abundance, density profile, and linear bias as functions of the linear power spectrum rms σ(r_0.2) and the MG free parameter (either f_R0 or z_SSB). We fit for the free coefficients of these parameterizations using the measurements made on sets of N-body simulations. Applying these parameterizations to analyze the same simulations, we recovered values for the MG free parameters that are compatible with the true values within 2σ in the case of a joint analysis including all void probes (abundance, density profile, and linear bias). Additionally, the values of the MG parameter coming independently from void abundance and void density profile analyses are also compatible with each other within two standard deviations.

The constraints on MG parameters from the linear bias are weaker than those from the density profile and abundance analyses, mainly because we analyzed relatively small box simulations, that is, cubic boxes with side 256 Mpc h⁻¹. This provides a poor sampling of Fourier modes on linear scales. Nevertheless, the constraints on MG parameters from the bias analysis show that we can distinguish between GR and a f(R) model with |f_R0| > 9.3 × 10⁻⁷ at the 95% confidence level. Similarly, we can distinguish between GR and a symmetron model with z_SSB >  0.8 at the 95% confidence level.

We also applied the symmetron analysis to the f(R) simulations and vice versa. For the MG scenarios closest to GR, that is, z_SSB = 1 and |f_R0| = 10⁻⁶, we were unable to significantly distinguish between symmetron and f(R) using any of the void properties we analyzed, even though we can distinguish them from GR. However, analyzing the simulations with an incorrect theory causes a difference in the MG parameter best fits inferred from individual probes, indicating that the MG model used is inappropriate.

For the other MG scenarios, z_SSB = {2,  3 } and |f_R0| = {10⁻⁵,  10⁻⁴}, we can distinguish among f(R), symmetron, and GR based mainly on two features. First, the MG parameters posterior distributions for the abundance and the density profile are compatible with each other within two standard deviations when the correct model is used, but with an incorrect theory, they are inconsistent by over four standard deviations. Second, the minimum χ²/d.o.f. is between 1.3 and 2.5 times larger when an incorrect theory is applied.

Finally, the joint analysis shows a difference of over three standard deviations between GR and the weakest modification on the MG models we analyzed. This type of analysis appears a promising tool for distinguishing gravity models, but further studies must be made, including realistic observational conditions. We expect the combination of void and halo properties to be particularly useful for constraining and distinguishing MG models. Because halos and voids respond differently to the increased forces and the screening effects that are unique to each model, the joint analysis of halo and void properties is expected to provide important consistency tests and help break degeneracies in parameter space. We hope to address some of these questions in future work.

Acknowledgments

We thank Claudio Llinares for providing the N-body simulations used in this work. This work has made use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by the Brazilian agency FAPESP (grant 2009/54006-4) and the INCT-A. EP and RV are supported by FAPESP. ML is partially supported by FAPESP and CNPq. DFM acknowledges support from the Research Council of Norway, and the NOTUR facilities.

References

Appendix A: Best fits

In this appendix we show the abundance, density profile, and bias associated with the best-fit MG parameters (see Figs. 4 and 5) recovered from the joint analysis of the three void properties.

Fig. A.1. Void density profiles measured in the f(R) simulations (points), and best fit of the model in Eq. (16) (lines). We split the void catalog into seven subsamples, corresponding to different void sizes, then we stacked all the void density-profiles for each subsample. Here we show the stacked subsamples with mean radius r0.2 = 3.0 Mpc h−1 (top panel), 5.2 Mpc h−1 (middle panel), and 9.0 Mpc h−1 (bottom panel).

Fig. A.2. Same as Fig. A.1, but for the symmetron simulations. In every case, the void wall is higher for stronger deviations of the modified gravity with respect to GR.

Fig. A.3. Linear matter-void bias for the f(R) and the symmetron cases. The points and error bars represent the measurements from simulations described in Sect. 4.3, while the solid lines represent the best fit of |fR0| or zSSB in the f(R) or symmetron cases, respectively.

Fig. A.4. Matter-void bias as a function of scale k for voids with different sizes given by the value of r0.2 in the ΛCDM simulation. Points denote simulation measurements and lines represent the best fit of the large-scale trend. A linear function in k2 was fit for large-scale modes with k <  0.25 Mpc h−1.

Fig. A.5. Similar to Fig. A.3, but for the void abundance measured and fit in f(R) and symmetron simulations, showing more large (fewer small) voids in the MG scenarios.

All Tables

All Figures

Fig. 1. Illustration in 2D of the assignment of void candidates to families (voids). We show voids that have already been identified {Vi,  Vj,  Vk,  Vl}, with both Vk and Vl having more than one member. The centers of each member in a family are inside the radius of another member belonging to a previous generation of the same family (we highlight this by showing the center of the family members, except for the very first detected member of each family). Blue circles denote void candidates and blue shades represent their cores (which correspond to 70% of the candidate’s radius).

In the text

	Fig. 2. Constraints on the f(R) parameter fR0 derived by the matter-void linear bias analysis.
In the text

	Fig. 3. Constraints on the symmetron parameter zSSB derived by the matter-void linear bias analysis.
In the text

	Fig. 4. Posterior distributions for the free parameter of the f(R) theory when we analyzed the {ΛCDM,  f6,  f5,  f4} simulations from top to bottom. We note that |fR0| constraints from abundance and density profile are compatible at {1,  2,  1,  1} σ for the {ΛCDM,  f6,  f5,  f4} cases; see Table 12.
In the text

	Fig. 5. Posterior distributions for the free parameter of the symmetron theory. The zSSB constraints from abundance and density profile are compatible at {1,  1,  2,  2} σ for the {ΛCDM,  A,  B,  D} cases, respectively, see Table 12.
In the text

	Fig. 6. Recovered value for |fR0| for the symmetron simulations.
In the text

	Fig. 7. Recovered value for zSSB for the f(R) simulations.
In the text

Fig. A.1. Void density profiles measured in the f(R) simulations (points), and best fit of the model in Eq. (16) (lines). We split the void catalog into seven subsamples, corresponding to different void sizes, then we stacked all the void density-profiles for each subsample. Here we show the stacked subsamples with mean radius r0.2 = 3.0 Mpc h−1 (top panel), 5.2 Mpc h−1 (middle panel), and 9.0 Mpc h−1 (bottom panel).

In the text

	Fig. A.2. Same as Fig. A.1, but for the symmetron simulations. In every case, the void wall is higher for stronger deviations of the modified gravity with respect to GR.
In the text

	Fig. A.3. Linear matter-void bias for the f(R) and the symmetron cases. The points and error bars represent the measurements from simulations described in Sect. 4.3, while the solid lines represent the best fit of |fR0| or zSSB in the f(R) or symmetron cases, respectively.
In the text

	Fig. A.4. Matter-void bias as a function of scale k for voids with different sizes given by the value of r0.2 in the ΛCDM simulation. Points denote simulation measurements and lines represent the best fit of the large-scale trend. A linear function in k2 was fit for large-scale modes with k <  0.25 Mpc h−1.
In the text

	Fig. A.5. Similar to Fig. A.3, but for the void abundance measured and fit in f(R) and symmetron simulations, showing more large (fewer small) voids in the MG scenarios.
In the text