Issue 
A&A
Volume 632, December 2019



Article Number  A52  
Number of page(s)  14  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201935949  
Published online  27 November 2019 
Cosmic voids in modified gravity scenarios
^{1}
Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern, 0315 Oslo, Norway
^{2}
Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, Rua do Matão 1371, São Paulo, SP 05508090, Brazil
email: elduartep@usp.br
Received:
24
May
2019
Accepted:
24
September
2019
Modified gravity (MG) theories aim to reproduce the observed acceleration of the Universe by reducing the dark sector while simultaneously recovering General Relativity (GR) within dense environments. Void studies appear to be a suitable scenario to search for imprints of alternative gravity models on cosmological scales. Voids cover an interesting range of density scales where screening mechanisms fade out, which reaches from a density contrast δ ≈ −1 close to their centers to δ ≈ 0 close to their boundaries. We present an analysis of the level of distinction between GR and two modified gravity theories, the Hu–Sawicki f(R) and the symmetron theory. This study relies on the abundance, linear bias, and density profile of voids detected in Nbody cosmological simulations. We define voids as connected regions made up of the union of spheres with a mean density given by ρ̅_{v} = 0.2 ρ̅_{m}, but disconnected from any other voids. We find that the height of void walls is considerably affected by the gravitational theory, such that it increases for stronger gravity modifications. Finally, we show that at the level of dark matter Nbody simulations, our constraints allow us to distinguish between GR and MG models with f_{R0} > 10^{−6} and z_{SSB} > 1. Differences of bestfit values for MG parameters that are derived independently from multiple void probes may indicate an incorrect MG model. This serves as an important consistency check.
Key words: largescale structure of Universe
© 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 socalled 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, longrange forces as well as cumulative effects due to different timeevolution paths can modify the largescale 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 intermediatesized 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 backgrounddense wall. Methods to detect them include freeparameter geometrical algorithms intended to find nonspherical 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 voidfinding 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 voidfinding algorithm, the void analysis can be performed in a selfconsistent 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 voidrelated 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 Nbody simulations we analyze in this work, as well as the voidfinder algorithm we employ. In Sect. 4 we describe the three voidrelated observables we choose to analyze, how we estimate them from the Nbody 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
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
and
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 ϕ:
The effective potential V_{eff}, Eq. (4), has a minimum at ϕ = 0 in highdensity environments, that is, where the local matter density ρ_{m} ≫ μ^{2}M^{2} ≡ ρ_{SSB}. In this case, V_{eff} holds the ϕ → −ϕ symmetry. On the other hand, in lowdensity 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):
which correspond to the range of the scalar field in Mpc h^{−1} units (λ_{0}), 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
where ℒ_{m} is the Lagrangian describing the matter fields, and
In this case, , where H_{0} 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_{1}/c_{2} = 6Ω_{Λ}/Ω_{m}.
This f(R) theory can be transformed into a scalar–tensor theory upon both the identification
and the frame transformation
In this case, the Compton wavelength or range of propagation of the scalar field at redshift zero is given by
where f_{R0} can be expressed as a function of {c_{2}, n} as
3. Methods
3.1. Nbody simulations
The cosmological Nbody 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^{3} dark matter particles in a (256 Mpc h^{−1})^{3} cubic box. The initial conditions correspond to a flat ΛCDM cosmology with parameters {Ω_{c}, Ω_{b}, h, σ_{8}, 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. Nbody 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^{−4}, 10^{−5}, 10^{−6}} for the f(R) theory, and by β = 1, λ_{0} = 1 Mpc h^{−1} 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.
Symmetron simulations analyzed in this work.
f(R) simulations analyzed in this work.
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. Voidfinding algorithm
Nonspherical 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 nonspherical voids also depends on the void size (Park & Lee 2007; Shoji & Lee 2012). Using the nonspherical version of our voidfinder 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 Einsteinde 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 voidfinding algorithm starts by computing the dark matter density on a regular grid of size l. This is done by applying the cloudincell (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^{−1} 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 voidcandidate 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 (dumbbellshaped) void, which is not the type of void we are interested in.
After all the spheres in the candidate catalog were cast, roughly onethird of them were discarded. The remaining twothirds of the original candidates were gathered into families (which we call nonspherical 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 familycasting 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 {V_{i}, V_{j}, V_{k}, V_{l}}, with both V_{k} and V_{l} 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. Nonspherical and spherical voids
The nonspherical 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 nonspherical 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 nonspherical members.
We extracted a sphericalvoid catalog from the nonspherical catalog by defining a spherical void as the largest member of each nonspherical 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 sphericalvoid catalog from a nonspherical one instead of directly detecting spherical voids from the very beginning weakens the problem of breaking highly nonspherical 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 voidmatter linear bias b(σ). Here, σ^{2} = σ^{2}(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 NavarroFrenkWhite 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 Nbody 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 voidmatter 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
which are slight modifications of linear functions in log_{10}f_{R0} or z_{SSB}. If they were pure linear functions (γ_{3} = γ_{6} = 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^{−8} instead of 0. On the other hand, γ_{4} could be written as γ_{4} = γ_{1} + γ_{2}log_{10}(γ_{3}) 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 1to1 map for f(R) and Symmetron is available, we were unable to find a common parameterization for the MG dependence. Furthermore, this nonequivalence 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),
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 voiddetection algorithm is not sensitive to the voidincloud process, which prevents small voids from surviving in dense regions (Jennings et al. 2013). From our point of view, the voidincloud 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 Nbody 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 voidincloud 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 PressSchechter 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 numbercounts 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 Nbody simulations. The errors on the measurements, estimated as the variance of the octantsubsamples 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
with
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:
where is the background dark matter density, δ_{0} 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:
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^{−1}, 14 Mpc h^{−1}) 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 − β0)α} with (1 − β_{0}) < 0 and α > 0, describing the tail of the profile. This tail is steeper for larger voids because it is the linear voidmatter 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 mattervoid 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 δ_{0} to a constant value because we did not take the inner part of the profile into account in this analysis. We computed δ_{0} 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. Mattervoid linear bias
The linear bias between the dark matter density field and the void distribution on large scales was estimated as
where P_{mm} is the matter power spectrum and P_{mv} is the mattervoid 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 shotnoise much higher than P_{vm}. Second, on large scales, P_{vm} changes sign for voids with r_{0.2} ∼ 7 Mpc h^{−1}, while P_{vv} is always positive by definition. This change of sign in b is expected by the peakbackgroundsplit (PBS) prediction (Chan et al. 2014) and has been seen in simulations (Hamaus et al. 2014b). Therefore, using P_{vm} to estimate the voidmatter 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 longwavelength background contribution, and a shortwavelength 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 selfconsistent linear halomatter 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 Nbody simulations including the MG effects.
We parameterize the linear voidmatter bias as
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^{−1}, 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^{−1} to fit the linear trend of the mattervoid 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^{2} for k < 0.25 Mpc h^{−1} (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 voidmatter 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.
Constraints on f_{R0} when the true value is f_{R0} = 0.
Constraints on f_{R0} when the true value is f_{R0} = 10^{−6}.
Constraints on f_{R0} when the true value is f_{R0} = 10^{−5}.
Constraints on f_{R0} when the true fiducial value is f_{R0} = 10^{−4}.
Constraints for f_{R0} when the fiducial value is z_{SSB} = 0.
Constraints for z_{SSB} when the fiducial value is z_{SSB} = 1.
Constraints for z_{SSB} when the fiducial value is z_{SSB} = 2.
Constraints for z_{SSB} when the fiducial value is z_{SSB} = 3.
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^{−6} 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 f_{R0} derived by the mattervoid linear bias analysis. 
Fig. 3. Constraints on the symmetron parameter z_{SSB} derived by the mattervoid 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 f_{R0} 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 z_{SSB} 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 highmass 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 voidmatter 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 tracermatter 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 χ^{2}/d.o.f. are higher when we assume an incorrect model.
Fig. 6. Recovered value for f_{R0} for the symmetron simulations. 
Fig. 7. Recovered value for z_{SSB} for the f(R) simulations. 
χ^{2} per degree of freedom for the joint analysis using all void observables.
indicates the level difference in the MG parameter recovered by the density and the abundance analyses with uncertainty σ.
6.1. Analyzing symmetron simulations using f(R) theory
In the case of a joint analysis using all void probes, the lowest χ^{2}/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 bestfit 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 χ^{2}/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 bestfit 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 χ^{2}/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^{−6}). The values of χ^{2}/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 socalled void selection function, described by completeness and purity functions, for example, which are specific to each voidfinder 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 lowmass 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^{−4} and z_{SSB} = 3), applying the correct MG theory analysis provides much better consistency between the different tests, and the minimum χ^{2}/d.o.f. is considerably lower (compared to the incorrect MG theory). Therefore at the level of Nboby simulations, we can safely distinguish between symmetron cases with z_{SSB} ≳ 2 and f(R) cases with f_{R0}≳10^{−5}.
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 sphericalvoid finding algorithm to construct void catalogs in Nbody 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 voidfinding 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. Nonspherical voids were defined as the union of spheres with a neighbortoneighbor overlap larger than a given threshold, while spherical voids were selected as the largest sphere of the nonspherical 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 Nbody 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^{−1}. 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^{−7} 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^{−6}, 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^{−5}, 10^{−4}}, 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 χ^{2}/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 Nbody 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/540064) and the INCTA. 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
 Adhikari, S., Sakstein, J., Jain, B., Dalal, N., & Li, B. 2018, JCAP, 1811, 033 [Google Scholar]
 Aghamousa, A., Aguilar, J., Ahlen, S., et al. 2016, ArXiv eprints [arXiv:1611.00036] [Google Scholar]
 Aguena, M., & Lima, M. 2018, Phys. Rev. D, 98, 123529 [NASA ADS] [CrossRef] [Google Scholar]
 Baker, T., Clampitt, J., Jain, B., & Trodden, M. 2018, Phys. Rev. D, 98, 023511 [NASA ADS] [CrossRef] [Google Scholar]
 Barreira, A., Cautun, M., Li, B., Baugh, C. M., & Pascoli, S. 2015, JCAP, 2015, 028 [CrossRef] [Google Scholar]
 Bolejko, K., Clarkson, C., Maartens, R., et al. 2013, Phys. Rev. Lett., 110, 021302 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Davis, A.C., Li, B., & Winther, H. A. 2012, Phys. Rev. D, 86, 044015 [NASA ADS] [CrossRef] [Google Scholar]
 Burrage, C., Copeland, E. J., Käding, C., & Millington, P. 2019, Phys. Rev. D, 99, 043539 [NASA ADS] [CrossRef] [Google Scholar]
 Cai, Y.C., Padilla, N., & Li, B. 2015, MNRAS, 451, 1036 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. C., Hamaus, N., & Desjacques, V. 2014, Phys. Rev. D, 90, 103521 [NASA ADS] [CrossRef] [Google Scholar]
 Colberg, J. M., Sheth, R. K., Diaferio, A., Gao, L., & Yoshida, N. 2005, MNRAS, 360, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Colberg, J. M., Pearce, F., Foster, C., et al. 2008, MNRAS, 387, 933 [NASA ADS] [CrossRef] [Google Scholar]
 Contigiani, O., Vardanyan, V., & Silvestri, A. 2019, Phys. Rev. D, 99, 064030 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Elyiv, A., Marulli, F., Pollina, G., et al. 2015, MNRAS, 448, 642 [NASA ADS] [CrossRef] [Google Scholar]
 Frenk, C. S., Cautun, M., & Cai, Y.C. 2016, MNRAS, 457, 2540 [NASA ADS] [CrossRef] [Google Scholar]
 Hagstotz, S., Gronke, M., Mota, D. F., & Baldi, M. 2019, A&A, 629, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamaus, N., Sutter, P. M., & Wandelt, B. D. 2014a, Phys. Rev. Lett., 112, 251302 [NASA ADS] [CrossRef] [Google Scholar]
 Hamaus, N., Wandelt, B. D., Sutter, P. M., Lavaux, G., & Warren, M. S. 2014b, Phys. Rev. Lett., 112, 041304 [NASA ADS] [CrossRef] [Google Scholar]
 Hamaus, N., Pollina, G., Weller, J., et al. 2017, MNRAS, 469, 787 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Hinterbichler, K., & Khoury, J. 2010, Phys. Rev. Lett., 104, 231301 [NASA ADS] [CrossRef] [Google Scholar]
 Hinterbichler, K., Khoury, J., Levy, A., & Matas, A. 2011, Phys. Rev. D, 84, 103521 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004 [CrossRef] [Google Scholar]
 Jennings, E., Li, Y., & Hu, W. 2013, MNRAS, 434, 2167 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, B. J. T., Van De Weygaert, R., & Platen, E. 2007, MNRAS, 380, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, A., Jain, B., Khoury, J., & Trodden, M. 2015, Phys. Rep., 568, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kowalski, M., Rubin, D., Aldering, G., et al. 2008, ApJ, 686, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K. 2016, Rep. Progr. Phys., 79, 046902 [CrossRef] [Google Scholar]
 Krause, E., Chang, T.C., Doré, O., & Umetsu, K. 2012, ApJ, 762, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Kreisch, C. D., Pisani, A., Carbone, C., et al. 2019, MNRAS, 488, 4413 [NASA ADS] [CrossRef] [Google Scholar]
 Lam, T. Y., Clampitt, J., Cai, Y.C., & Li, B. 2015, MNRAS, 450, 3319 [NASA ADS] [CrossRef] [Google Scholar]
 Lavaux, G., & Wandelt, B. D. 2012, ApJ, 754, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Zhao, G., & Koyama, K. 2012, MNRAS, 421, 3481 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., & Mota, D. F. 2014, Phys. Rev. D, 89, 084023 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., Mota, D. F., & Winther, H. A. 2014, A&A, 562, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lopes, R. C. C., Voivodic, R., Abramo, L. R., & Sodré, Jr., L. 2018, JCAP, 1809, 010 [NASA ADS] [CrossRef] [Google Scholar]
 Lopes, R. C., Voivodic, R., Abramo, L. R., & Sodré, Jr., L. 2019, JCAP, 2019, 026 [CrossRef] [Google Scholar]
 Maggiore, M., & Riotto, A. 2010, ApJ, 711, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Nadathur, S., Hotchkiss, S., Diego, J. M., et al. 2014, Proc. Int. Astron. Union, 11, 542 [CrossRef] [Google Scholar]
 Nadathur, S., Carter, P. M., Percival, W. J., Winther, H. A., & Bautista, J. E. 2019, Phys. Rev. D, 100, 023504 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Neyrinck, M. C. 2008, MNRAS, 386, 2101 [NASA ADS] [CrossRef] [Google Scholar]
 Olive, K. A., & Pospelov, M. 2008, Phys. Rev. D, 77, 043524 [NASA ADS] [CrossRef] [Google Scholar]
 Park, D., & Lee, J. 2007, Phys. Rev. Lett., 98, 081301 [NASA ADS] [CrossRef] [Google Scholar]
 Perlmutter, S. 2003, Phys. Today, 56, 53 [CrossRef] [Google Scholar]
 Perlmutter, S., Aldering, G., Valle, M. D., et al. 1998, Nature, 391, 51 [CrossRef] [Google Scholar]
 Planck Collaboration VI. 2019, A&A, in press [arXiv:1807.06209] [Google Scholar]
 Pogosian, L., & Silvestri, A. 2010, Phys. Rev. D, 81, 049901 [NASA ADS] [CrossRef] [Google Scholar]
 Rapetti, D. A., Morris, R. G., Allen, S. W., et al. 2008, MNRAS, 383, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Ricciardelli, E., Quilis, V., & Planelles, S. 2013, MNRAS, 434, 1192 [NASA ADS] [CrossRef] [Google Scholar]
 Ricciardelli, E., Quilis, V., & Varela, J. 2014, MNRAS, 440, 601 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Sánchez, C., Clampitt, J., Kovacs, A., et al. 2017, MNRAS, 465, 746 [NASA ADS] [CrossRef] [Google Scholar]
 Scaramella, R., Mellier, Y., Amiaux, J., et al. 2014, IAU Symp., 306, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F., Lima, M., Oyaizu, H., & Hu, W. 2009, Phys. Rev. D, 79, 083518 [NASA ADS] [CrossRef] [Google Scholar]
 Schuster, N., Hamaus, N., Pisani, A., et al. 2019, Arxiv eprints [arXiv:1905.00436] [Google Scholar]
 Sheth, R. K., & Van De Weygaert, R. 2004, MNRAS, 350, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Shoji, M., & Lee, J. 2012, ArXiv eprints [arXiv:1203.0869] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Tsujikawa, S. 2008, Phys. Rev. D, 77, 023507 [NASA ADS] [CrossRef] [Google Scholar]
 van de Weygaert, R., & van Kampen, E. 1993, MNRAS, 263, 481 [NASA ADS] [Google Scholar]
 Voivodic, R., Lima, M., Llinares, C., & Mota, D. F. 2017, Phys. Rev. D, 95, 024018 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., & Ferreira, P. G. 2015, Phys. Rev. D, 92, 064005 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., Mota, D. F., & Li, B. 2012, ApJ, 756, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., Schmidt, F., Barreira, A., et al. 2015, MNRAS, 454, 4208 [NASA ADS] [CrossRef] [Google Scholar]
 Wyman, M., Jennings, E., & Lima, M. 2013, Phys. Rev. D, 88, 084029 [NASA ADS] [CrossRef] [Google Scholar]
 Yoo, J., & Watanabe, Y. 2012, Int. J. Mod. Phys. D, 21, 1230002 [NASA ADS] [CrossRef] [Google Scholar]
 Zivick, P., Sutter, P. M., Wandelt, B. D., Li, B., & Lam, T. Y. 2015, MNRAS, 451, 4215 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Best fits
In this appendix we show the abundance, density profile, and bias associated with the bestfit 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 densityprofiles for each subsample. Here we show the stacked subsamples with mean radius r_{0.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 mattervoid 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 f_{R0} or z_{SSB} in the f(R) or symmetron cases, respectively. 
Fig. A.4. Mattervoid bias as a function of scale k for voids with different sizes given by the value of r_{0.2} in the ΛCDM simulation. Points denote simulation measurements and lines represent the best fit of the largescale trend. A linear function in k^{2} was fit for largescale 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
indicates the level difference in the MG parameter recovered by the density and the abundance analyses with uncertainty σ.
All Figures
Fig. 1. Illustration in 2D of the assignment of void candidates to families (voids). We show voids that have already been identified {V_{i}, V_{j}, V_{k}, V_{l}}, with both V_{k} and V_{l} 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 f_{R0} derived by the mattervoid linear bias analysis. 

In the text 
Fig. 3. Constraints on the symmetron parameter z_{SSB} derived by the mattervoid 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 f_{R0} 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 z_{SSB} 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 f_{R0} for the symmetron simulations. 

In the text 
Fig. 7. Recovered value for z_{SSB} 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 densityprofiles for each subsample. Here we show the stacked subsamples with mean radius r_{0.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 mattervoid 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 f_{R0} or z_{SSB} in the f(R) or symmetron cases, respectively. 

In the text 
Fig. A.4. Mattervoid bias as a function of scale k for voids with different sizes given by the value of r_{0.2} in the ΛCDM simulation. Points denote simulation measurements and lines represent the best fit of the largescale trend. A linear function in k^{2} was fit for largescale 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.