Cosmic voids in modified gravity scenarios

TL;DR

This paper investigates how cosmic voids can be used to differentiate between General Relativity and modified gravity theories by analyzing void properties in cosmological simulations, revealing potential observational signatures.

Contribution

It introduces a novel analysis of void abundance, bias, and density profiles to distinguish between GR and MG models, focusing on the Hu-Sawicki $f(R)$ and symmetron theories.

Findings

Void wall height increases with stronger gravity modifications.

Constraints can distinguish MG models with |f_{R0}| > 10^{-6}.

Void properties serve as consistency checks for MG theories.

Abstract

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 $\delta \approx -1$ close to their centers to $\delta \approx 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 n-body cosmological simulations. We define voids as connected regions made up of the union of spheres with a {\it \textup{mean}} density given by $\overline\rho_v=0.2\,\overline\rho_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 n-body simulations, our constraints allow us to distinguish between GR and MG models with $|f_{R0}| > 10^{-6}$ and $z_{SSB} > 1$. Differences of best-fit 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.

Tables

Table 1. Table 1: Symmetron simulations analyzed in this work. All cases have λ 0 = 1 subscript 𝜆 0 1 \lambda_{0}=1 Mpc/ h ℎ h and β = 1 𝛽 1 \beta=1 .

Symmetron case	$\Lambda$CDM	A	B	D
$z_{SSB}$	0	1	2	3

Table 2. Table 2: f ​ ( R ) 𝑓 𝑅 f(R) simulations analyzed in this work. All cases have n = 1 𝑛 1 n=1 and β = 1 / 6 𝛽 1 6 \beta=1/\sqrt{6} , and λ 0 = λ 0 ​ ( | f R ​ 0 | ) subscript 𝜆 0 subscript 𝜆 0 subscript 𝑓 𝑅 0 \lambda_{0}=\lambda_{0}(|f_{R0}|) is given by Eq. ( 10 ) in Mpc/ h ℎ h units.

$f\left(R\right)$ case	$\Lambda$CDM	f6	f5	f4
$\left|f_{R0}\right|$	0	${10}^{-6}$	${10}^{-5}$	${10}^{-4}$
${\lambda }_0\left(\left|f_{R0}\right|\right)$	0	2.4	7.5	23.7

Table 3. Table 3: Constraints on | f R ​ 0 | subscript 𝑓 𝑅 0 |f_{R0}| when the true value is | f R ​ 0 | = 0 subscript 𝑓 𝑅 0 0 |f_{R0}|=0 . For the GR simulation the f ​ ( R ) 𝑓 𝑅 f(R) analysis shows that log 10 ⁡ | f R ​ 0 | < − 6.03 subscript 10 subscript 𝑓 𝑅 0 6.03 \log_{10}|f_{R0}|<-6.03 at the 95 % percent 95 95\% of confidence level when the bias analysis is applied.

${\log }_{10}\left|f_{R0}\right|$	best fit	mean	1$\sigma$	2$\sigma$	3$\sigma$
abundance	-8.54	-7.88	¡-7.71	¡-7.36	¡-7.14
density profile	-10.00	-8.68	¡-8.54	¡-8.17	¡-7.96
bias	-10.00	-6.80	¡-6.55	¡-6.03	¡-5.66
all	-10.00	-8.72	¡-8.57	¡-8.20	¡-7.98

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\frac{R}{2}{M}_{pl}^2-\frac{{\nabla }_{\mu }\varphi {\nabla }^{\mu }\varphi }{2}-V(\varphi )\right]+S_m({\tilde{g}}_{\mu \nu },{\psi }_i),$$

$$A(\varphi )=1+\frac{1}{2}{\left(\frac{\varphi }{M}\right)}^2,$$

$$V(\varphi )=V_0-\frac{{\mu }^2{\varphi }^2}{2}+\frac{\lambda {\varphi }^4}{4}.$$

$$\square \varphi =V_{\text{eff}}=\frac{1}{2}\left(\frac{{\rho }_m}{{\mu }^2}-M^2\right){\varphi }^2+\frac{\lambda {\varphi }^4}{4}.$$

$${\lambda }_0=\frac{1}{\sqrt{2}\mu },\beta =\frac{{\varphi }_0{M}_{pl}}{M^2},(1+z_{SSB}{)}^3=\frac{{\mu }^2{M}^2}{{\overline{\rho }}_{m0}},$$

$$S=\int d^{4}x\sqrt{-\tilde{g}}\left[\frac{\tilde{R}+f(\tilde{R})}{16\pi G}+{\mathcal{L}}_m\right],$$

$$f(\tilde{R})\equiv -\frac{m^2{c}_1(\tilde{R}/m^2{)}^n}{c_2(\tilde{R}/m^2{)}^n+1}.$$

$$f_R\equiv \frac{df(\tilde{R})}{d\tilde{R}}=e^{-2\beta \varphi /M_{pl}}-1\approx -\frac{2\beta \varphi }{M_{pl}},$$

$${\tilde{g}}_{\mu \nu }=e^{2\beta \varphi /M_{pl}}{g}_{\mu \nu },\text{with}\beta =\frac{1}{\sqrt{6}}.$$

$${\lambda }_0=3\sqrt{\frac{n+1}{{\Omega }_m+4{\Omega }_{\Lambda }}}\sqrt{\frac{|f_{R0}|}{10^{-6}}}\frac{\text{Mpc}}{h},$$

$$f_{R0}\equiv {f_R|}_{z=0}=-\frac{6n{\Omega }_{\Lambda }}{c_2{\Omega }_m}{\left(\frac{{\Omega }_m/3}{{\Omega }_m+4{\Omega }_{\Lambda }}\right)}^{n+1}.$$

$$\tilde{\gamma}(x)=\left\{\begin{array}[]{ll}\gamma_{1}+\gamma_{2}\,\log_{10}(\gamma_{3}+x)\,,&x\equiv|f_{R0}|\text{ for $f(R)$}\,,\\ \gamma_{4}+\gamma_{5}\,\cfrac{x}{1+\gamma_{6}\,e^{-x^{2}}}\,,&x\equiv z_{SSB}\text{ for symmetron}\,,\end{array}\right.$$

$$\begin{array}{rl}\frac{dn}{d\ln R}(\sigma ,x)=& \frac{f(\sigma ,x)}{V(R)}\frac{d\ln {\sigma }^{-1}}{d\ln R},\\ f(x,\sigma )=& A{\sigma }^{\tilde{\gamma }(x)}(1+{\nu }^b)e^{-c{\sigma }^{-2}},\nu =1.686/\sigma ,\end{array}$$

$${\rho }_v(r)=\frac{3m}{4\pi }\sum _i\frac{\Theta (r_i;r,\delta r)}{(r+\delta r{)}^3-(r-\delta r{)}^3},$$

$$\Theta(r_{i};r,\delta r)=\left\{\begin{array}[]{ll}1\,,&\text{ when }r_{i}\in(r-\delta r,r+\delta r)\,,\\ 0\,,&\text{ otherwise }\,,\end{array}\right.$$

$$\frac{{\rho }_v(r)}{{\overline{\rho }}_m}-1={\delta }_0\frac{1-G(yc{)}^{\alpha }}{1+(yc{)}^{\beta }},y=r/r_{0.2},$$

$$\begin{array}{rl}\alpha (x,\sigma )=& \tilde{\alpha }(x)-{\alpha }_0\sigma ,\\ \beta (x,\sigma )=& {\beta }_0\alpha (x,\sigma ),\\ c^{\beta (x,\sigma )}(\sigma )=& c_0{\sigma }^{-c_1},\\ G(x,\sigma )=& \tilde{G}(x)+G_0{\log }_{10}(\sigma ),\end{array}$$

$$b={\frac{P_{mv}(k)}{P_{mm}(k)}|}_{k\to 0},$$

$$b(r_{0.2})=a+\tilde{c}(x){\sigma }^{-2}+d{\sigma }^{-4},\sigma =\sigma (r_{0.2}),$$

Full text

††offprints: E. L. D. Perico,

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 $\Lambda$CDM (Perlmutter et al., 2003; Eisenstein et al., 2005; Kowalski et al., 2008; Rapetti et al., 2008; Hinshaw et al., 2013; Planck Collaboration et al., 2018). 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 $\Lambda$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 $\Lambda$, 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; Sánchez et al., 2017; Nadathur et al., 2019), and gravitational lensing (Sánchez et al., 2016; 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 $\Lambda$CDM outcome.

Specifically, we aim to distinguish between the $\Lambda$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 section 2 we give an overview of the MG models we are interested in. In section 3 we describe the n-body simulations we analyze in this work, as well as the void-finder algorithm we employ. In section 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 sections 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 section 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 $\phi$. The action for a scalar-tensor theory can be written as

[TABLE]

where $S_{m}$ is the action associated with the standard matter fields $\psi_{i}$. The fields $\psi_{i}$ are coupled to the scalar field $\phi$ through the Jordan frame metric $\tilde{g}_{\mu\nu}$, which is related to the Einstein frame metric $g_{\mu\nu}$ by $\tilde{g}_{\mu\nu}=A^{2}(\phi)g_{\mu\nu}$. For the particular case of the symmetron theory, we have

[TABLE]

and

[TABLE]

In this case, both $A(\phi)$ and $V(\phi)$ are symmetric upon the transformation $\phi\rightarrow-\phi$. Here, $\mu$ and $M$ are mass scales, $\lambda$ is a dimensionless coupling constant, $M_{pl}^{2}=(8\pi G)^{-1}$ is the Planck mass scale, and $g$ is the determinant of $g_{\mu\nu}$.

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

[TABLE]

The effective potential $V_{\text{eff}}$, Eq. (4), has a minimum at $\phi=0$ in high-density environments, that is, where the local matter density $\rho_{m}\gg\mu^{2}M^{2}\equiv\rho_{SSB}$. In this case, $V_{\text{eff}}$ holds the $\phi\rightarrow-\phi$ symmetry. On the other hand, in low-density environments, $\rho_{m}\ll\rho_{SSB}$, the effective potential displays two minima at $\phi=\pm\phi_{0}\sqrt{1-\rho_{m}/\rho_{SSB}}$, breaking the symmetry of the model. Here $\phi_{0}=\mu/\sqrt{\lambda}$ is the expected value of $\phi$ for $\rho_{m}=0$.

The free parameters of the symmetron model $\{\mu,\,M,\,\lambda\}$ are usually exchanged by the physical parameters associated with the scalar field for $\rho_{m}=0$ (Winther et al., 2012):

[TABLE]

which correspond to the range of the scalar field in Mpc$/h$ units ($\lambda_{0}$), the dimensionless coupling strength to matter ($\beta$), 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 $\rho_{SSB}$, for which the symmetry breaking takes place through the expression $\rho_{SSB}=\overline{\rho}_{m0}(1+z_{SSB})^{3}$. Here $\overline{\rho}_{m0}$ 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

[TABLE]

where $\mathcal{L}_{m}$ is the Lagrangian describing the matter fields, and

[TABLE]

In this case, $m^{2}=H_{0}^{2}\Omega_{m0}$, where $H_{0}$ and $\Omega_{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\Omega_{\Lambda}/\Omega_{m}$.

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

[TABLE]

and the frame transformation

[TABLE]

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

[TABLE]

where $f_{R0}$ can be expressed as a function of $\{c_{2},\,n\}$ as

[TABLE]

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, R., 2002) to include MG models, with $512^{3}$ dark matter particles in a $(256\,\text{Mpc}/h)^{3}$ cubic box. The initial conditions correspond to a flat $\Lambda$CDM cosmology with parameters $\{\,\Omega_{c}\,,\,\Omega_{b}\,,\,h\,,\,\sigma_{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, Steffen 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^{-4},\,10^{-5},\,10^{-6}\}$ for the $f(R)$ theory, and by $\beta=1$, $\lambda_{0}=1\leavevmode\nobreak\ \text{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 $\Lambda$CDM or GR case, where no gravity modifications are included. We assume here that the $\Lambda$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 $\Delta_{\text{v}}=0.2$ that reaches as far as the void radius (see Appendix A of Jennings et al. (2013)). Even though $\Delta_{\text{v}}$ depends on cosmology and MG parameters, we set it to its EdS value as it 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 $\Delta_{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 $\rho(r_{j})\equiv 0.2\times\overline{\rho}_{m}$, which is the expected density for spherical voids at redshift $z=0$ (see Jennings et al., 2013). Here, $\overline{\rho}_{m}$ 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 $\@vec{x}_{j}$ (initially the center of the cell), but also at the corners of a cube of side $d$ around $\@vec{x}_{j}$. ii) When one of the corners maximizes the size of the sphere that is refined, we moved $\@vec{x}_{j}$ to that position. Otherwise we did not update $\@vec{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 $\@vec{x}_{j}$ of a given sphere $j$ were refined, we set $r_{k}=0$ for every grid sphere $k$ whose center $\@vec{x}_{k}$ was closer than $0.9r_{k}$ to $\@vec{x}_{j}$. This was done to avoid duplicates because under refinement, these grid spheres $k$ would converge to the same values $r_{j}$ and $\@vec{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 cold dark matter (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}$.

3.3 Non-spherical and spherical voids

The non-spherical void catalog consists of all the families that were founded by applying steps 1) to 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 $0.2\times\overline{\rho}_{m}$ 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(\sigma)$, the void density profile $\rho(r\,,r_{0.2})$, and the void-matter linear bias $b(\sigma)$. Here, $\sigma^{2}=\sigma^{2}(R)$ is the variance of the linear matter power spectrum smoothed on a scale $R$, which was computed for the $\Lambda$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 $\sigma$ 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

[TABLE]

which are slight modifications of linear functions in $\log_{10}|f_{R0}|$ or $z_{SSB}$. If they were pure linear functions ($\gamma_{3}=\gamma_{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 $\sim 10^{-8}$ instead of [math]. On the other hand, $\gamma_{4}$ could be written as $\gamma_{4}=\gamma_{1}+\gamma_{2}\log_{10}(\gamma_{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 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 $\delta_{\text{v}}$ and $\delta_{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),

[TABLE]

where $\{\,A\,,\,b\,,\,c\,,\,\tilde{\gamma}\,\}$ 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 $\tilde{\gamma}(x)$, 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. 12 we show the measured void abundance and the best fit of Eq. (13) for the $\Lambda$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

[TABLE]

with

[TABLE]

where $m$ and $r_{i}$ are the mass and position of the dark matter particles, while $2\delta 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 $\rho(r)\to\rho(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:

[TABLE]

where $\bar{\rho}_{m}$ is the background dark matter density, $\delta_{0}$ is a constant, and the parameters $\{\,\alpha\,,\,\beta\,,\,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 $\sigma(r_{0.2})$ as follows:

[TABLE]

$\tilde{\alpha}(x)$ and $\tilde{G}(x)$ have the same functional form as $\tilde{\gamma}(x)$ 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), $\sigma\equiv\sigma(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}\in(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).

$\alpha\,,\beta$) For $r\gtrsim 2\,r_{0.2}$ the expression (16) falls like $r^{\alpha-\beta}$, in our case as $r^{(1-\beta_{0})\alpha}$ with $(1-\beta_{0})<0$ and $\alpha>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 $\sigma$ dependence of $\alpha$. The MG dependence in the term $\tilde{\alpha}(x)$ 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^{-\beta}$ 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^{-\beta}$ increase with $\sigma$ (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, $\tilde{G}(x)$ 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. 8 and 9. We show $\rho_{v}(r/r_{0.2})/\bar{\rho}_{m}-1$ for three different void sizes and the different MG cases.

We set the density at the void center $\delta_{0}$ to a constant value because we did not take the inner part of the profile into account in this analysis. We computed $\delta_{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 Matter-void linear bias

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

[TABLE]

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 $b=\sqrt{P_{vv}/P_{mm}}$. 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}\sim 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

[TABLE]

where $a$ and $d$ are constants, while $\tilde{d}=\tilde{c}(x)$ has the same functional form as $\tilde{\gamma}(x)$ 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. 11 for the $\Lambda$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$ $h$/Mpc (see Fig. 11).

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. Figs. 4 and 5 and Tables 5-5 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}|\gtrsim 10^{-6}$ or $z_{SSB}\gtrsim 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.

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. 12.

• •

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. 9. 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 section 6.