# Chameleon Screening in Cosmic Voids

Andrius Tamosiunas  Chad Briddon  Clare Burrage  Alan Cutforth  Adam Moss  and
Thomas Vincent

###### Abstract

A key goal in cosmology in the upcoming decade will be to form a better
understanding of the accelerated expansion of the Universe. Upcoming surveys,
such as the Vera C. Rubin Observatory’s 10-year Legacy Survey of Space and
Time (LSST), Euclid and the Square Killometer Array (SKA) will deliver key
datasets required to tackle this and other puzzles in contemporary cosmology.
With this data, constraints of unprecedented power will be put on different
models of dark energy and modified gravity. In this context it is crucial to
understand how screening mechanisms, which hide the deviations of these
theories from the predictions of general relativity in local experiments,
affect structure formation. In this work we approach this problem by using a
combination of analytic and numerical methods to describe chameleon screening
in the context of cosmic voids. We apply a finite element method (FEM) code,
SELCIE, to solve the chameleon equation of motion for a number of void
profiles derived from observational data and simulations. The obtained results
indicate a complex relationship between the properties of cosmic voids and the
size of the chameleon acceleration of a test particle. We find that the fifth
force on a test particle in a void is primarily related to the depth and the
inner density gradient of the void. For realistic void profiles, the obtained
chameleon-to-Newtonian acceleration ratios range between $a_{\phi}/a_{\rm
Newt}\approx 10^{-6}-10^{-5}$. However, it should be noted that in unusually
deep voids with large inner density gradients, the acceleration ratios can be
significantly higher. Similarly, other chameleon models, such as $f(R)$ Hu-
Sawicki theory allow for significantly higher acceleration ratios. Given these
results, we also discuss the optimal density profiles for detecting the fifth
force in the upcoming observational surveys.

##  1 Introduction

We are living in the era of the accelerated expansion of the Universe, as
illustrated by the wealth of data from type Ia supernovae, baryon acoustic
oscillations (BAO), galaxy clusters and other observational probes [1, 2, 3,
4]. Arguably the simplest and the most successful explanation for the
accelerating expansion as of yet is the cosmological constant. Nonetheless,
the discrepancy between the value of the cosmological constant determined in
cosmology and that calculated in quantum field theory could hint towards a
more complicated model for the accelerated expansion. Similarly, contemporary
problems in cosmology, such as the $H_{0}$ tension, could be hinting towards
new physics. In this context, modifying general relativity offers an
alternative explanation to some of the puzzles outlined above [5, 6, 7, 8, 9].

A subset of theories of special interest are scalar field models that possess
a screening mechanism. These are theories that turn off modified gravity
effects in and around high density regions, hence satisfying the tight Solar
System observational constraints, while still possibly allowing interesting
effects on cosmological scales. Well known examples of screening mechanisms
include the chameleon, the symmetron and the Vainshtein models [10, 11, 12,
13]. Despite the elusive nature of screened fifth forces mediated by these
types of scalar fields, there exists a number of techniques to detect them in
laboratory and cosmology tests. In particular, extensive cosmological and
laboratory constrains have been put on chameleon screening [14, 15, 16, 17,
18]. Similarly, techniques for constraining symmetron and Vainshtein screening
are well described in the literature as well [19, 20, 21, 22, 23, 24].

In the context of cosmological tests of gravity, a system of special interest
is that of cosmic voids. Voids are the largest underdensities in the Universe
and their properties are inherently linked to the physics of large scale
structure. Given these properties, cosmic voids make ideal cosmological
laboratories for testing screened gravity theories. Using voids as a testing
ground for models of modified gravity has been explored extensively in the
literature (e.g. see Ref. [25, 26, 27, 28]). Nonetheless, many unanswered
questions regarding the properties of cosmic voids remain to be addressed.
Specifically, an interesting question to ask is how does the fifth force
depend on the density distribution in cosmic voids? Previous work has been
done to deduce the magnitude of the fifth force in voids; however, these
analyses have been performed primarily using a top hat/step function density
models for the voids (e.g. see Ref. [29]). Recent simulations and
observational studies indicate that void density distributions are
significantly more complicated than that of a step function. Hence, one might
ask what is the relation between different void density distributions and the
magnitude of the fifth force? As an extension of our previous numerical work
(Ref. [18]), here we attempt to answer this question by numerically solving
the chameleon equation of motion (EOM) for different void density profiles.

In the following sections we introduce the model at hand along with the
numerical methods used to solve the EOM. Specifically, section 2 introduces
the specific chameleon model that we will investigate. Section 3 briefly
describes the SELCIE code used to solve the chameleon EOM. Section 4
introduces some basic properties of voids, with a particular emphasis on
different density profiles used in our analysis. The analytic and the
numerical results are described in sections 5 and 6. Appendix A shows a
comparison between the results obtained in this work and the corresponding
results for certain $f(R)$ theories. In this work we use natural units
($c=\hbar=1$) and a mostly positive metric signature
$\eta_{\mu\nu}=\mathrm{diag}(-1,1,1,1)$. We denote the reduced Planck mass as
$M_{\rm Pl}=(8\pi G)^{-1}$, with $G$ as the Newton’s constant. In all the
shown figures we used rescaled radial units, such that $\hat{r}=r/L$, with
$L=2\times 31.68$ Mpc/$h$ (twice the typical void size as described in Ref.
[30]). Similarly, the densities are rescaled by the critical density of the
Universe $\hat{\rho}=\rho/\rho_{\mathrm{c}}$ with
$\rho_{\mathrm{c}}=3H^{2}/8\pi G=1.8788\times
10^{-26}h^{2}\mathrm{~{}kg}\mathrm{~{}m}^{-3}$ with $h=0.673$.

##  2 Chameleon Gravity

A non-minimally coupled scalar-tensor theory is described by the following
action:

$S=\int\mathrm{d}x^{4}\sqrt{-g}\left(\frac{M_{\mathrm{Pl}}^{2}}{2}R-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi-V(\phi)\right)-\int\mathrm{d}x^{4}\mathcal{L}_{m}\left(\varphi_{m}^{(i)},\tilde{g}_{\mu\nu}^{(i)}\right),$
(2.1)

with $\phi$ as the scalar field, $V(\phi)$ as the potential, $\mathcal{L}_{m}$
as the matter Lagrangian, $\varphi_{m}^{(i)}$ as the matter fields and
$\tilde{g}_{\mu\nu}^{(i)}$ as the Jordan frame metric. The superscript i here
refers to the i-th matter species. The Jordan frame metric can be related to
the Einstein frame metric via the following relation:

$\tilde{g}_{\mu\nu}^{(i)}=A_{i}^{2}(\phi)g_{\mu\nu}.$ (2.2)

Chameleon gravity refers to a particular choice of $A_{i}$ and $V(\phi)$, e.g.

$A_{i}(\phi)=e^{\beta_{i}\phi/M_{Pl}},$ (2.3)

with $\beta_{i}$ term here as the coupling constant for the i-th matter
species and

$V(\phi)=\Lambda^{4}\left(1+\frac{\Lambda^{n}}{\phi^{n}}\right),$ (2.4)

with an energy scale $\Lambda$ (note that other choices of the potential that
lead to chameleon screening are possible). Assuming a universal coupling to
matter, the EOM is given by:

$\nabla^{2}\phi=-\frac{n\Lambda^{n+4}}{\phi^{n+1}}+\frac{\rho}{M}.$ (2.5)

Note that here we have also assumed a time-independent scalar field and
defined $M\equiv M_{\mathrm{Pl}}/\beta$, where we dropped the index for
$\beta$ under the assumption of a universal coupling constant.

Figure 1: Observational and laboratory constraints of chameleon gravity. The
regions excluded by each specific test are indicated in the figure. The figure
on the left shows the constraints on $\Lambda$ and $M$, for $n=1$. The figure
on the right shows the corresponding constraints for $n$ and $M$ for
$\Lambda=2.4$ meV. The region labelled astrophysics contains the bounds from
both Cepheid and rotation curve tests. The dashed line indicates the dark
energy scale $\Lambda=2.4$ meV. The black, red, and blue arrows in both
figures show the lower bound on $M$ coming from neutron bouncing and
interferometry experiments. In this work we focus on the white regions of the
parameter space. Figure adapted from [16].

For the above potential, Eq. 2.4, the different chameleon models correspond to
different choices of the three parameters: $n$, $\Lambda$, $M$. In this work
we primarily focus on the models that are still allowed by the current
observational constraints (see Fig. 1). Specifically, we investigate the
parameter space regions still allowed by the observational constraints to
determine which region would produce the highest values for the chameleon
acceleration. In accordance with our previous work in Ref. [18], we will find
that in order to maximise the chameleon acceleration, the smallest allowed
value of $M$ is needed. Similarly, the smallest allowed value of $n$ is
required to maximise the chameleon acceleration. Lastly, we want $\Lambda$ to
be as large as possible. Putting all these requirements together the parameter
values that maximise the fifth force, while still being in the allowed
parameter space, are: $n=1$, $M=10^{-14}\times M_{\rm Pl}$ and
$\Lambda=10^{-4}$ eV.

Here we should also note that there are other chameleon models often studied
in the literature, such as the $f(R)$ Hu-Sawicki model (e.g. Ref. [31]), which
corresponds to a different region of the parameter space shown in Fig 1.
Studying such chameleon models in detail is left for future work, however, an
illustrative comparison between our results and the corresponding results in
the Hu-Sawicki model is given in appendix A.

The chameleon EOM is non-linear and finding exact analytic solutions for
complex non-symmetric density distributions is generally not possible. Hence,
to find solutions we will apply the numerical technique known as FEM. FEM
allows finding solutions to non-linear equations for arbitrary density
distributions. Specifically, in order to solve the EOM, we will employ the
software package SELCIE, which is described further in the upcoming section.

##  3 The SELCIE Solver

FEM refers to a widely used numerical technique for solving linear and non-
linear differential equations in one, two or three spatial dimensions. FEM has
now existed for over 80 years and its range of applications varies widely from
fluid mechanics and engineering to meteorology and physics [32, 33, 34].
Recently FEM has also been applied extensively in the context of modified
gravity [35, 36, 37, 18].

FEM subdivides the problem domain into smaller subdomains referred to as
finite elements. Specifically, the domain is segmented into cells, whose
boundaries are defined by their vertices. The value of the field inside each
cell is then approximated by a piecewise polynomial function that matches the
field values at each of the cell’s vertices. For non-linear second order
differential equations one can use Green’s theorem and Taylor expand the non-
linear term. The static solution can then be obtained at each mesh point using
algorithms such as the Picard or the Newton iteration methods.

In the case of the chameleon EOM, we start by rescaling the equation
(following our previous approach in Ref. [37]):

$\alpha\hat{\nabla}^{2}\hat{\phi}=-\hat{\phi}^{-(n+1)}+\hat{\rho},$ (3.1)

where $\hat{\phi}=\phi/\phi_{\infty}$, $\hat{\rho}=\rho/\rho_{\infty}$,
$\hat{\nabla}^{2}=L^{2}\nabla^{2}$, $L$ is the domain size and all the
constants have been absorbed into the $\alpha$ parameter. The $\rho_{\infty}$
and $\phi_{\infty}$ here refer to the background density and the corresponding
background value of the field. The $\alpha$ parameter is given by:

$\alpha\equiv\left(\frac{M\Lambda}{L^{2}\rho_{\infty}}\right)\left(\frac{nM\Lambda^{3}}{\rho_{\infty}}\right)^{\frac{1}{n+1}}.$
(3.2)

To linearise Eq. 3.1 we can Taylor expand the non-linear term around some
estimate of the field, $\phi_{k}$:

$\begin{gathered}\hat{\phi}^{-(n+1)}\approx\hat{\phi}_{k}^{-(n+1)}-(n+1)\hat{\phi}_{k}^{-(n+2)}\left(\hat{\phi}-\hat{\phi}_{k}\right)+\mathcal{O}\left(\hat{\phi}-\hat{\phi}_{k}\right)^{2}\\\
\approx(n+2)\hat{\phi}_{k}^{-(n+1)}-(n+1)\hat{\phi}_{k}^{-(n+2)}\hat{\phi}+\mathcal{O}\left(\hat{\phi}-\hat{\phi}_{k}\right)^{2}.\end{gathered}$
(3.3)

We can then obtain the required variational (integral) form for the EOM, Eq.
3.1, by applying Green’s theorem and multiplying each term by $v_{j}$, which
is an arbitrary test function that vanishes on the boundary of our domain,
$\partial\Omega$. Finally, substituting the expansion for the non-linear term,
Eq. 3.3, into the variational form of the EOM, we obtain:

$\alpha\int_{\Omega}\hat{\nabla}\hat{\phi}\cdot\hat{\nabla}v_{j}\mathrm{~{}d}x+\int_{\Omega}(n+1)\hat{\phi}_{k}^{-(n+2)}\hat{\phi}v_{j}\mathrm{~{}d}x=\int_{\Omega}(n+2)\hat{\phi}_{k}^{-(n+1)}v_{j}\mathrm{~{}d}x-\int_{\Omega}\hat{\rho}v_{j}\mathrm{~{}d}x.$
(3.4)

We can solve the non-linear EOM iterativly, by solving the above linearised
version, Eq. 3.4, for some estimate of the field $\hat{\phi}_{k}$ and updating
the estimate using

$\hat{\phi}_{k+1}=\omega\hat{\phi}+(1-\omega)\hat{\phi}_{k},$ (3.5)

until the solution has converged to the desired degree of accuracy. The
$\omega$ in Eq. 3.5 is the relaxation parameter, which controls the step size
between the consecutive solutions found by the solver. This iterative
procedure is called the Picard method.

To automate the solution to Eq. 3.4, we employ the SELCIE algorithm111SELCIE
is available at: https://github.com/C-Briddon/SELCIE. [37]. SELCIE is a FEM
algorithm that automates the solution to the chameleon EOM for arbitrary
density distributions. It allows for an easy mesh creation procedure as well
as using different solvers to obtain solutions in 2D and 3D. SELCIE is based
on the FEniCS Project software package, which is a collection of free and
open-source software modules dedicated to automating solutions to differential
equations via the FEM [38, 39, 40, 41]. Note that the approach laid out here
is nearly identical to the methods previously applied to study chameleon
gravity in NFW halos [18]. The SELCIE software package and the procedure of
solving the chameleon EOM is described in detail in Ref. [37, 18]. In summary,
to obtain the solutions, we used the FEniCS in-built preconditioned Krylov
solver with the conjugate gradient solver method and the default
preconditioner settings (see section 4.3 in Ref. [37] for further details).
The mesh precision (the number of mesh cells per unit length) was set to 150.
The code execution was set to stop when the absolute change in the field
variable becomes smaller than $\delta\hat{\phi}=10^{-14}$. The relaxation
parameter, $\omega$, was set to 0.3, a value determined by experimentation.

##  4 Cosmic Voids

Cosmic voids are the most underdense regions in the Universe, found between
the network of filaments in the cosmic web. The astrophysical properties of
voids make them objects of special interest in cosmology and studies of
gravity. Namely, being a major component of the cosmic web, as well as some of
the largest objects in the Universe, voids hold a wealth of information about
the underlying cosmology. The voids’ structures, shapes and mutual alignments
are sensitive to the properties of dark energy and gravity. Voids, being
nearly empty (yet still populated by galaxies), also form a natural
environment to study galaxy formation and evolution. Finally, the low
densities found in voids make them a great environment for studying screening
mechanisms.

The properties of cosmic voids have been studied extensively, both
observationally and by using large scale structure formation simulations.
Estimates of void sizes range between 10-50 $\mathrm{Mpc}/h$ [42, 43]. A more
recent study based on the Cosmic Void Catalog discovered similar results, with
most void sizes falling between 10-30 $\mathrm{Mpc}/h$ [44]. While the vast
majority of the voids have sizes near the middle of the mentioned range, there
are examples of extremely large/small voids. As an example, the Bootes void is
know as one of the largest objects in the Universe with an approximate radial
size of 44-50 $\mathrm{Mpc}/h$ [45, 46]. On the lower side of the size
spectrum, mini-voids have been detected with sizes in the range of 0.7-3.5
$\mathrm{Mpc}/h$ [47, 48].

There have been a number of studies investigating the density profiles of
cosmic voids. Multiple universal void profiles have been proposed, based on
both simulations and observational data. The simplest type of a profile
explored in the literature is the step function (top hat) profile, described
by:

$\hat{\rho}_{\rm ST}(\hat{r})=\begin{cases}\hat{\rho}_{\rm
out},&\hat{r}>\hat{r}_{\rm step}\\\ \hat{\rho}_{\rm
in},&\hat{r}\leq\hat{r}_{\rm step}\end{cases}.$ (4.1)

Here $\hat{\rho}_{\rm in}$ and $\hat{\rho}_{\rm out}$ refer to the inner and
the outer density and $\hat{r}_{\rm step}$ is the radial position of the step.
Note that the densities here are rescaled by the critical density of the
Universe, while the radial variables are divided by the size of the domain of
the problem at hand, $L$. While such a profile is simple to work with, and has
been explored in the literature extensively, e.g. [29], in this work we will
show that such a profile is not realistic enough to model the chameleon models
that are left unconstrained by the current observational data. A more
realistic profile choice is based on the hyperbolic tangent function, given
by:

$\hat{\rho}_{\rm th}(\hat{r})=\frac{1}{2}\bigg{[}(\hat{\rho}_{\rm
in}+\hat{\rho}_{\rm out})+(\hat{\rho}_{\rm out}-\hat{\rho}_{\rm
in})\tanh(k(\hat{r}-\hat{r}_{\rm step}))\bigg{]}.$ (4.2)

This profile converges to a Heaviside step function for $k\rightarrow\infty$.
As before, the densities and the radial variables are rescaled.

In Ref. [49] a more physically motivated four-parameter profile based on
observational data and simulation studies is introduced:

$\hat{\rho}_{\mathrm{v}}(\hat{r})-1=\delta_{\mathrm{c}}\frac{1-\left(\hat{r}/\hat{r}_{\mathrm{s}}\right)^{\alpha_{\mathrm{v}}}}{1+\left(\hat{r}/\hat{r}_{\mathrm{v}}\right)^{\beta_{\mathrm{v}}}}+\gamma_{\mathrm{v}}.$
(4.3)

The $\delta_{\mathrm{c}}$ parameter here is the central density contrast,
$\hat{r}_{\mathrm{v}}$ is the void radius rescaled by the domain size,
$\hat{r}_{\mathrm{s}}$ is the rescaled scale radius at which
$\rho_{\mathrm{v}}=\bar{\rho}$ (or $\hat{\rho}_{v}=1$ in rescaled units) while
$\alpha_{\mathrm{v}}$ and $\beta_{\mathrm{v}}$ determine the inner and outer
slope of the void density distribution. The extra parameter
$\gamma_{\mathrm{v}}$ is introduced in some works to account for the fact that
a given void can be located in a under/over-dense region (e.g. as discussed in
Ref. [30]).

In Ref. [50, 51] a variation of the profile in Eq. 4.3 is fitted to simulated
void data. Specifically, the void radius is replaced with the scale radius in
the denominator of Eq. 4.3 and fitted to realistic mock luminous red galaxy
(LRG) catalogues from the Jubilee simulation, as well as void catalogues
constructed from the SDSS LRG and Main Galaxy samples [52, 53, 51]. The voids
were detected using the watershed transform algorithm ZOBOV [54]. The self-
similarity and universality of void density profiles were further studied in
Ref. [55], while the velocity profiles were studied in Ref. [56]. In Ref. [30]
the profile described in Eq. 4.3 is fitted to the lensing potential map from
Planck CMB lensing reconstruction analysis data in combination with the
“Public Cosmic Void Catalog”. The analysis in Ref. [30] determined that
introducing an extra parameter, $\gamma_{\mathrm{v}}$, improves the fit to the
data. This is done to account for the fact that Sloan Digital Sky Survey voids
reside in an underdense region. In our analysis we will employ the profiles
and the best-fit numbers from Ref. [50, 30] as a benchmark to compare the
calculated results against. Specifically, for Ref. [50], the following best-
fit numbers will be used: $\delta_{\mathrm{c}}=-0.69$,
$\hat{r}_{\mathrm{s}}=0.81\hat{r}_{\mathrm{v}}$, $\alpha_{\mathrm{v}}=1.57$,
$\beta_{\mathrm{v}}=5.72$. Similarly, for the best-fit profile from Ref. [30],
we will use $\delta_{\mathrm{c}}=-0.78$,
$\hat{r}_{\mathrm{s}}=1.26\hat{r}_{\mathrm{v}}$, $\alpha_{\mathrm{v}}=2.69$,
$\beta_{\mathrm{v}}=13.85$. For the void size we will use a value of
$r_{\mathrm{v}}=31.68$ Mpc/$h$, which is a good estimate for the size of a
typical void as discussed in Ref [30]. The $\gamma_{\mathrm{v}}$ parameter
will generally be set to zero, but, for the sake of completeness, the effects
of the mentioned parameter will be investigated in Fig. 7. The referenced
profiles are plotted in Fig. 2, where we see that there is significant
variation between the best-fit void profiles in different works in the
literature. This is not entirely surprising, as the best-fit parameters in
Ref. [50, 30] have been determined using different techniques and come with
different systematics and errors. In this work this will not be an issue, as
we will only use the outlined profiles as a benchmark, and the void parameters
will be varied to explore a large set of possible profiles.

Figure 2: Void density profiles determined from the CMB and SDSS data and
simulations as discussed in [30] and [50]. The density is rescaled by the
value of the critical density of the Universe, while the radial value is
divided by the size of the domain: $L=2\times 31.68\;\mathrm{Mpc}/h$, i.e.
twice the typical size of a void. The values in the brackets refer to the
best-fit parameters: ($\delta_{\rm c}$, $\alpha_{\rm v}$, $\beta_{\rm v}$,
$\hat{r}_{\rm s}/\hat{r}_{\rm v}$, $\gamma_{\rm v}$). The blue band
corresponds to the standard deviation around the best-fit values.

Given the form of Eq. 4.3, one can derive an analytical expression for the
void mass enclosed by some radius $R$ assuming spherical symmetry:

$\begin{split}M_{\mathrm{v}}(R)=\int_{0}^{R}4\pi
r^{2}\rho_{\mathrm{v}}(r)dr&=4\pi\rho_{c}R^{3}\Bigg{[}\frac{1}{3}+\frac{\delta_{c}}{3(\alpha_{\mathrm{v}}+3)}\bigg{(}(\alpha_{\mathrm{v}}+3)\;{}_{2}F_{1}\bigg{(}1,\frac{3}{\beta_{\mathrm{v}}};\frac{\beta_{\mathrm{v}}+3}{\beta_{\mathrm{v}}};-\bigg{(}\frac{R}{r_{\mathrm{v}}}\bigg{)}^{\beta_{\mathrm{v}}}\bigg{)}\\\
&\quad-3\bigg{(}\frac{R}{r_{\mathrm{s}}}\bigg{)}^{\alpha_{\mathrm{v}}}{}_{2}F_{1}\bigg{(}1,\frac{\alpha_{\mathrm{v}}+3}{\beta_{\mathrm{v}}};\frac{\alpha_{\mathrm{v}}+\beta_{\mathrm{v}}+3}{\beta_{\mathrm{v}}};\bigg{(}\frac{R}{r_{\mathrm{s}}}\bigg{)}^{\beta_{\mathrm{v}}}\bigg{)}\bigg{)}+\frac{\gamma_{\mathrm{v}}}{3}\Bigg{]},\end{split}$
(4.4)

where ${}_{2}F_{1}$ corresponds to the hypergeometric function, defined as:

${}_{2}F_{1}(a,b;c;z)=\sum_{n=0}^{\infty}\frac{(a)_{n}(b)_{n}}{(c)_{n}}\frac{z^{n}}{n!}=1+\frac{ab}{c}\frac{z}{1!}+\frac{a(a+1)b(b+1)}{c(c+1)}\frac{z^{2}}{2!}+\cdots$
(4.5)

Integrating the Poisson equation then gives the relation for the regular
Newtonian gravitational acceleration:

$a_{\mathrm{Newt}}(R)=-\frac{GM_{\mathrm{v}}(R)}{R^{2}}.$ (4.6)

##  5 Analytic Investigation

When it comes to analyzing the effects of cosmic void density profiles on
chameleon screening, a lot of headway can be made analytically. As previously
discussed in Ref. [18], for much of the allowed parameter space (when dealing
with cosmological length scales) the chameleon field tracks the minimum value
of its effective potential. More specifically, in terms of the EOM (Eq. 2.5),
this refers to the regime where $\alpha\nabla^{2}\phi\ll V^{\prime}_{\rm
eff}(\phi)$, with $V_{\rm eff}(\phi)=V(\phi)+\rho\phi/M$. We will refer to
this as the small-$\alpha$ regime. In this regime we can use the
approximation:

$\hat{\phi}(\hat{r})\approx[\hat{\rho}_{\mathrm{v}}(\hat{r})]^{-\frac{1}{n+1}}.$
(5.1)

An approximate analytic expression for the field gradient is therefore given
by the derivative of Eq. 5.1:

$\frac{\partial\hat{\phi}(\hat{r})}{\partial\hat{r}}\approx-\frac{[\hat{\rho}_{\rm
v}(\hat{r})]^{-\frac{n+2}{n+1}}}{n+1}\frac{\partial\hat{\rho}_{\rm
v}(\hat{r})}{\partial\hat{r}}.$ (5.2)

Eq. 5.2 illustrates that the chameleon field depends strongly on the shape of
the void density profile. This behaviour is illustrated in Fig. 3. In
particular, we plot different density profiles, starting with the universal
density profile described in Fig. 2 and ending with a step function (top hat)
profile. For the more realistic Hamaus-Chantavat-type profile (Eq. 4.3), we
also vary the $\alpha_{\mathrm{v}}$ parameter (which controls the shape of the
inner slope of the void), while keeping the other parameters fixed to their
best-fit values. A general conclusion that can be made here is that the more
step function-like a profile is, the larger the corresponding radial
derivative of the chameleon field. The radial derivative of the field is
maximised for the step function profile in the region of the step. In this
case, the radial derivative of the density and the corresponding chameleon
field goes to infinity. For this reason, the analytic approximation for the
small-$\alpha$ regime, Eq. 5.1, cannot be used to accurately describe the
behaviour of the field at the step region. The full numerical solutions that
we will describe in Section 6 are required to compute the chameleon field at
the step region.

While the step function profile is clearly not a realistic description of the
density distribution inside of a void, it is an interesting limiting case. To
investigate this limiting case further we will explore the hyperbolic tangent
profile described in Eq. 4.2. Such a profile converges to a step function for
large values of $k$. The resulting ratios between the chameleon and the
Newtonian acceleration are shown in Fig. 4. Specifically, the acceleration
ratios are shown for different values of $k$ along with the different depths
of the void. The key finding here is that the acceleration ratios are
primarily limited by the inner slope of the void profile along with the inner
depth of the density distribution. In principle one could obtain arbitrarily
large values for the acceleration ratio by increasing the inner slope of the
density distribution. As before, however, for large values of $k$ the
small-$\alpha$ approximation breaks down and full numerical treatment is
needed. By applying a small perturbation to the field value given by Eq. 5.1,
we can find that the small-$\alpha$ approximation will hold only if
$\sqrt{\alpha}\ll 1/k$. For $\alpha$ to lie in the allowed regions of the
parameter space it must satisfy $\alpha<10^{-12}$. Therefore, so long as $k\ll
10^{6}$, the small-$\alpha$ approximation will hold for all values of $\alpha$
allowed by the parameter space.

Figure 3: The results of the analytic investigation. Left: different density
profiles. The black dashed and the blue solid lines correspond to the profiles
described in Fig. 2. The dashed and the dotted blue lines were generated by
varying the $\alpha_{\rm v}$ parameter value, while keeping the other
parameters fixed to the best-fit values used for the solid blue line (i.e. the
Ref. [30] values). The dotted line corresponds to the step function density
profile. The arrow corresponds to the Dirac delta function (i.e. the value
goes off to infinity). Centre: the corresponding radial derivative of the
different density profiles. Right: the corresponding radial gradients of the
chameleon field. Note that all quantities are rescaled as defined in section
3.

Figure 4: The results of the analytic investigation for a tanh density
profile. Left: different density profiles for different values of the $k$
parameter. Center: the corresponding chameleon-to-Newtonian acceleration
ratios. Right: acceleration ratios for different inner void densities (i.e.
different depths). In this case each profile is calculated with $k=40$ and
$\hat{\rho}_{\mathrm{out}}=1.0$. In all cases the dotted black line
corresponds to the step function profile, while the black arrow denotes the
profile going infinite (Dirac delta function). Note that all quantities are
rescaled as defined in section 3.

##  6 Numerical Investigation

Figure 5: Void density profiles for different free parameters. The
$\alpha_{\mathrm{v}}$, $\beta_{\mathrm{v}}$ and
$\hat{r}_{s}/\hat{r}_{\mathrm{v}}$ parameters are varied, while keeping the
other parameters fixed to their best-fit values in each case (i.e. we are
varying the best-fit parameter values of the blue profile in Fig. 2). As
before the black dotted line and the blue band refers to the best-fit profiles
from Ref. [30] and [50]. Note that all quantities are rescaled as defined in
section 3.

Figure 6: The chameleon-to-Newtonian acceleration ratios for density profiles
depicted in Fig. 5.

In this section we turn our attention to more realistic density distributions.
In particular, we explore the Nadathur and Hamaus-Chantavat-type density
profiles from Ref. [50, 30] (as depicted in Fig. 2). Each free parameter is
varied (while keeping the other parameters fixed to their best-fit values).
For each set of the parameters we solve the EOM numerically by employing
SELCIE, as described in more detail in section 3. Fig. 5 depicts the various
density profiles generated by varying the $\alpha_{\mathrm{v}}$,
$\beta_{\mathrm{v}}$ and the $\hat{r}_{\mathrm{s}}/\hat{r}_{\mathrm{v}}$
parameters, while Fig. 6 shows the corresponding chameleon-to-Newtonian
acceleration ratios. Fig. 5 and 8 show the corresponding results for voids of
varying depths.

We find that for the more realistic profiles, as depicted in Fig. 5, the
acceleration ratio values are generally lower than in the case of the step
function and the tanh density profiles, Eq. 4.2. A key finding is that the
parameters that control the depth of the void along with the inner density
slope (i.e. $\delta_{\mathrm{c}}$ and $\alpha_{\mathrm{v}}$) have the most
significant effect on the values of the acceleration ratios. Similarly, the
$\gamma_{\mathrm{v}}$ parameter has a similar effect to the
$\delta_{\mathrm{c}}$ parameter through its effect on the central depth of the
void.

An interesting observation that can be drawn from the $\alpha_{\mathrm{v}}$
variation plot (Fig. 6) is that there are generally two regimes of behaviour.
Specifically, for $\alpha_{\mathrm{v}}\lesssim 2$ the acceleration ratios
reach maximal values towards the centre of the void. However, for
$\alpha_{\mathrm{v}}\gtrsim 2$, the acceleration ratio curve falls to some
maximally negative value before turning back towards zero in the central
region of the void. Such behaviour is determined by the gradient of the
density profile in the central region of the void along with the slope closer
to the region where the void density reaches its maximum value, i.e. the
compensation wall. On the other hand, varying $\beta_{\mathrm{v}}$ and the
$\hat{r}_{\mathrm{s}}/\hat{r}_{\mathrm{v}}$ ratio only changes the slope in
the outer regions of the void, thus the two regimes of behaviour are not
observed in those cases. The depth of the void is also of key importance for
determining the maximum value of the acceleration ratio. Generally, the deeper
the void, the larger the chameleon-to-Newtonian acceleration ratio. This can
be achieved by either having a small value of $\delta_{\mathrm{c}}$ or a large
negative value of $\gamma_{\mathrm{v}}$.

Figure 7: Void density profiles for different $\delta_{\rm c}$ (void depths)
and $\gamma_{\rm v}$ parameters. As before, the other parameters are set to
the best-fit values from [30] and [50]. Note that all quantities are rescaled
as defined in section 3.

Figure 8: Chameleon-to-Newtonian acceleration ratios corresponding to the
density profiles depicted in Fig. 7.

##  7 Conclusions

In this work we have investigated the surviving chameleon models in the
context of cosmic voids. Specifically, we have studied chameleon models still
allowed by the observational and laboratory constraints and solved the
chameleon EOM given a selection of void density profiles using both analytic
and numerical techniques. In the case of a step function profile, we showed
that the chameleon acceleration vanishes in all regions except for the step
itself. At the step, the radial derivative of the density profile goes to
infinity, which, under the small-$\alpha$ approximation, would result in an
infinite chameleon acceleration. However, as discussed in section 5, the
small-$\alpha$ approximation is not valid in profiles of sufficient steepness,
hence full numerical treatment is required, which results in large (but
finite) values at the step region. In summary, while the step-function profile
is not physically realistic, it is nonetheless useful as it acts as a limiting
case that maximises the chameleon acceleration. In other words, the more step
function-like a void is, the higher is the corresponding chameleon
acceleration. While, in principle, the void acceleration is only limited by
the steepness of the inner density gradient, in practice, for voids of
realistic gradients, chameleon-to-Newtonian acceleration ratios of
$a_{\phi}/a_{\rm Newt}\sim 10^{-6}-10^{-5}$ can be obtained for the parameter
values studied in this work.

When investigating the more physically motivated profiles we have shown that
the chameleon field profiles exhibit a rich and complex behaviour. In summary,
the chameleon-to-Newtonian acceleration ratios are most sensitive to the inner
density slope and the depth of the void. The variation of the inner density
slope, in particular, can lead to two different behaviours – the acceleration
ratios reaching maximum (negative) values towards the centre of the void or
somewhere between the center and the compensation wall of the void. However,
irrespective of which of these two regimes we consider, the acceleration
ratios reach similar values of $\sim 10^{-6}$.

It should be noted that chameleon models are not limited to the scalar-tensor
theories studied in this work. In particular, a subset of $f(R)$ theories,
such as Hu-Sawicki gravity, also exhibit chameleon screening [31]. It can be
shown that $f(R)$ theories are equivalent to scalar-tensor theories with a
particular choice for the potential (e.g. see Ref. [57]). The Hu-Sawicki
theory, in particular, corresponds to different scalar-tensor parameter
values, which lead to acceleration ratios that are much higher. The
approximate acceleration ratios for the Hu-Sawicki model are calculated in
appendix A. Currently with our code we can only calculate approximate results
for the Hu-Sawicki model. This is the case, as the code cannot handle the
negative fractional values for the $n$ parameter. To obtain a direct
comparison, a new implementation of the code would have to be written, with a
specific FEM approach applied to solve the $f(R)$ EOM directly. This is left
for future work.

Our results shine light on the observational prospects of detecting the
effects of the chameleon field in cosmic voids. For the parameter space region
studied in this work, the acceleration ratios of $\approx 10^{-6}-10^{-5}$
would be extremely challenging to detect. Nonetheless, even such a small fifth
force could potentially be detected statistically through its effects on
structure formation. A fifth force, in the central region of a void, or near
the compensation wall, could potentially increase the average void depth and
the size of the typical compensation wall. To investigate such effects in
detail, large scale structure simulations with modified gravity effects would
be required. Ideally, one would run hydrodynamic modified gravity simulations
to investigate the possible degeneracies between astrophysical effects and
effects due to a fifth force. Other models, such as the mentioned Hu-Sawicki
model are known to exhibit much larger fifth forces, which have been and will
be further constrained by previous and upcoming observational data (e.g.
galaxy cluster weak lensing measurements) [15, 58].

The study described in this work can be extended in multiple ways. For
instance, here we had made a key assumption of voids being spherically
symmetric. Real cosmic voids are better described as triaxial objects. While
generally voids tend to become more spherical with expansion, their axis
ratios vary between $\sim 0.5-1.0$ [48]. Hence, analogously to NFW halos,
which we had previously studied in Ref. [18], we expect fifth force effects
that depend on the void shape. Nonetheless, we do not expect these effects to
significantly increase or reduce the chameleon-to-Newtonian acceleration
ratios. Effects due to void shapes could also be investigated in detail by
employing modified gravity simulations.

Finally, it should be noted that in this work we had only investigated simple
void systems – i.e. single voids of different depths. There are, however,
other more complicated systems that could be investigated numerically. As an
example, a system of multiple voids separated by filaments would lead to a
large gradient in the density distribution, which, in turn, would lead to
relatively high values of the chameleon accleration. Such systems could be
investigated with SELCIE; however, a more detailed study would require
hydrodynamical modified gravity simulations. On a similar note, such systems
are also of special interest in other modified gravity models, e.g. the
symmetron, which are currently less constrained. The symmetron-type models
exhibit complex time-dependent behaviour (e.g. domain walls). Such behaviour
could be investigated by employing the SELCIE code; however, at the time of
writing, the code does not allow for the solving of time-dependent equations.
We are currently working towards expanding the code functionality to include a
solver for the symmetron EOM and to model time-dependent systems.

##  Appendix A Acceleration Ratios for Gravitational Strength Couplings

The acceleration ratios obtained in this work are significantly smaller that
those found in the literature (e.g. Ref. [59, 60]), where the fifth force is
found to be of similar magnitude to the Newtonian force inside cosmic voids.
The key reason for this discrepancy in results is due to the fact that the
mentioned works consider a different model, i.e. the $f(R)$ Hu-Sawicki model
[31]. The Hu-Sawicki model can be mapped directly to the parameter space of
the scalar-tensor chameleon models explored in this work. Namely, the Hu-
Sawicki model corresponds to a different part of the $\\{n,M,\Lambda\\}$
parameter space, resulting in significantly higher fifth force. In the rest of
this section, we lay out the mapping between the scalar-tensor models
described in this work and the corresponding $f(R)$ theories. Furthermore, we
calculate the acceleration ratios for a particular case of $f(R)$ gravity,
which is comparable to the Hu-Sawicki model.

The Hu-Sawicki model refers to a specific case of $f(R)$ gravity with:

$f(R)=-m^{2}\frac{c_{1}\left(R/m^{2}\right)^{n_{1}}}{c_{2}\left(R/m^{2}\right)^{n_{1}}+1},$
(A.1)

where $m,n_{1},c_{1}$ and $c_{2}$ are the model parameters. Our scalar-tensor
model can be directly compared to $f(R)$ by noting that:

$f_{R}=-\sqrt{\frac{2}{3}}\frac{\phi}{M_{\mathrm{Pl}}},$ (A.2)

with $f_{R}\equiv\partial f/\partial R$ and $\phi$ as the chameleon scalar
field. For the current value of the $f_{R}$ parameter we have:

$f_{R0}=-\sqrt{\frac{2}{3}}\frac{\phi_{\infty}}{M_{\mathrm{Pl}}},$ (A.3)

with $\phi_{\infty}$ as the field value at $r\rightarrow\infty$. Explicitly,
the background value of the field is:

$\phi_{\infty}=\left(\frac{nM\Lambda^{n+4}}{\rho_{\infty}}\right)^{\frac{1}{n+1}},$
(A.4)

where $\rho_{\infty}$ is the corresponding density. Substituting the field
value into Eq. A.3 gives:

$\left|f_{R0}\right|=\sqrt{\frac{2}{3}}\frac{1}{M_{\mathrm{Pl}}}\left[\frac{nM\Lambda^{n+4}}{\rho_{\infty}}\right]^{\frac{1}{n+1}}.$
(A.5)

Now we should note that the coupling parameter in $f(R)$ gravity is set to
$\beta=\sqrt{1/6}$ or equivalently $M=\sqrt{6}M_{\mathrm{Pl}}$. Furthermore,
we should note that models with $|f_{R0}|=\\{10^{-4},10^{-5},10^{-6}\\}$ are
often considered in the literature. Finally, if we set $n=1$222Here we should
note that the choice of $n=1$ does not directly correspond to the Hu-Sawicki
model. For a direct comparison we would need to set $n$ to a negative
fractional value, which currently our code cannot handle. Nonetheless, after
experimenting with different values, we found that the results do not strongly
depend on $n$. Hence our results should be a good approximation of the
acceleration ratios in the Hu-Sawicki model. for convenience and
$\rho_{\infty}=\rho_{c}$ we can then solve for $\Lambda$, i.e. the $\Lambda$
value that corresponds to the outlined $f_{R0}$ values in the corresponding
$f(R)$ model:

$\Lambda=\frac{3\rho_{\mathrm{c}}M_{\mathrm{Pl}}\left|f_{R0}\right|^{2}}{\sqrt{24}}.$
(A.6)

Figure 9: The acceleration ratios for $f(R)$ gravity with $|f_{R0}|=10^{-5}$.
The ratios were calculated using the best-fit parameter values for the density
profiles from Ref. [30] as shown in Fig. 2, with the $\alpha_{\mathrm{v}}$ and
$\delta_{\mathrm{c}}$ parameters varied correspondingly.

Numerically, if we plug in $\rho_{\rm c}=8.07\times 10^{-11}h^{2}$
$\mathrm{eV^{4}}$, $|f_{R0}|=10^{-5}$ and $M_{\rm Pl}=2.435\times 10^{27}$ eV,
we get $\Lambda=22.6$ eV. Finally, if we plug in the $n$, $M$ and $\Lambda$
values outlined above to our numerical code, it allows us to compute the
approximate values for the acceleration ratios in Hu-Sawicki gravity. For this
we used the Chantavat et al. (2017) best-fit numbers (Fig. 2) and varied the
$\alpha_{\rm v}$ and $\delta_{\rm c}$ parameters. This results in the
acceleration ratios given in Fig. 9. The obtained values are significantly
higher than those reported in the previous figures. In particular, the
acceleration ratios generally agree with the previous results in the
literature (e.g. Ref. [60]). Specifically, Table II in Ref. [60] lists the
acceleration ratio of $0.19\pm 0.002$ for void sizes in the range of
$25-35\;\mathrm{Mpc}/h$, redshift $z=0.5$ and $|f_{R0}|=10^{-5}$. Meanwhile,
our results vary between $|a_{\phi}/a_{\rm Newt}|\approx 0.1-1.0$ depending on
the density profile used. Here it should be noted that minor differences
between the results are to be expected as the acceleration values were
calculated using different datasets, corresponding to voids at different
redshifts, of different sizes and different density profiles.

## Acknowledgments

We would like to thank Seshadri Nadathur for his insights into different void
density profiles and various void finding algorithms.

Clare Burrage and Andrius Tamosiunas are supported by a Research Leadership
Award from The Leverhulme Trust. Chad Briddon is supported by the University
of Nottingham. Adam Moss is supported by a Royal Society University Research
Fellowship.

## References

  * [1] A.G. Riess, A.V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P.M. Garnavich et al., _Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant_ , _AJ_ 116 (1998) 1009 [astro-ph/9805201].
  * [2] V. Sahni and A. Starobinsky, _The Case for a Positive Cosmological $\Lambda$-Term_, _International Journal of Modern Physics D_ 9 (2000) 373 [astro-ph/9904398].
  * [3] P.J. Peebles and B. Ratra, _The cosmological constant and dark energy_ , _Reviews of Modern Physics_ 75 (2003) 559 [astro-ph/0207347].
  * [4] P. Ruiz-Lapuente, _Dark Energy: Observational and Theoretical Approaches_ , Cambridge University Press (2010).
  * [5] A. Joyce, L. Lombriser and F. Schmidt, _Dark Energy Versus Modified Gravity_ , _Annual Review of Nuclear and Particle Science_ 66 (2016) 95 [1601.06133].
  * [6] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, _Modified gravity and cosmology_ , _Phys. Rep._ 513 (2012) 1 [1106.2476].
  * [7] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri et al., _In the realm of the Hubble tension-a review of solutions_ , _Classical and Quantum Gravity_ 38 (2021) 153001 [2103.01183].
  * [8] M. Braglia, M. Ballardini, F. Finelli and K. Koyama, _Early modified gravity in light of the H 0 tension and LSS data_, _Phys. Rev. D_ 103 (2021) 043528 [2011.12934].
  * [9] D. Benisty and A.-C. Davis, _Dark energy interactions near the Galactic Center_ , _Phys. Rev. D_ 105 (2022) 024052 [2108.06286].
  * [10] A. Vainshtein, _To the problem of nonvanishing gravitation mass_ , _Phys. Lett. B_ 39 (1972) 393.
  * [11] J. Khoury and A. Weltman, _Chameleon cosmology_ , _Phys. Rev. D_ 69 (2004) 044026 [astro-ph/0309411].
  * [12] K. Hinterbichler, J. Khoury, A. Levy and A. Matas, _Symmetron cosmology_ , _Phys. Rev. D_ 84 (2011) 103521 [1107.2112].
  * [13] E. Babichev and C. Deffayet, _An introduction to the Vainshtein mechanism_ , _Classical and Quantum Gravity_ 30 (2013) 184001 [1304.7240].
  * [14] A. Terukina, L. Lombriser, K. Yamamoto, D. Bacon, K. Koyama and R.C. Nichol, _Testing chameleon gravity with the Coma cluster_ , _J. Cosmology Astropart. Phys_ 2014 (2014) 013 [1312.5083].
  * [15] H. Wilcox, D. Bacon, R.C. Nichol, P.J. Rooney, A. Terukina, A.K. Romer et al., _The XMM Cluster Survey: testing chameleon gravity using the profiles of clusters_ , _MNRAS_ 452 (2015) 1171 [1504.03937].
  * [16] C. Burrage and J. Sakstein, _Tests of chameleon gravity_ , _Living Reviews in Relativity_ 21 (2018) 1 [1709.09071].
  * [17] D.O. Sabulsky, I. Dutta, E.A. Hinds, B. Elder, C. Burrage and E.J. Copeland, _Experiment to detect dark energy forces using atom interferometry_ , _Phys. Rev. Lett._ 123 (2019) 061102.
  * [18] A. Tamosiunas, C. Briddon, C. Burrage, W. Cui and A. Moss, _Chameleon Screening Depends on the Shape and Structure of NFW Halos_ , _arXiv e-prints_ (2021) arXiv:2108.10364 [2108.10364].
  * [19] C. Llinares and L. Pogosian, _Domain walls coupled to matter: The symmetron example_ , _Phys. Rev. D_ 90 (2014) 124041 [1410.2857].
  * [20] C. Burrage, A. Kuribayashi-Coleman, J. Stevenson and B. Thrussell, _Constraining symmetron fields with atom interferometry_ , _J. Cosmology Astropart. Phys_ 2016 (2016) 041 [1609.09275].
  * [21] J. Beltrán Jiménez, F. Piazza and H. Velten, _Evading the Vainshtein Mechanism with Anomalous Gravitational Wave Speed: Constraints on Modified Gravity from Binary Pulsars_ , _Phys. Rev. Lett._ 116 (2016) 061101 [1507.05047].
  * [22] A. Dima and F. Vernizzi, _Vainshtein screening in scalar-tensor theories before and after GW170817: Constraints on theories beyond Horndeski_ , _Phys. Rev. D_ 97 (2018) 101302 [1712.04731].
  * [23] A. Hammami and D.F. Mota, _Probing modified gravity via the mass-temperature relation of galaxy clusters_ , _A &A_ 598 (2017) A132 [1603.08662].
  * [24] C. Llinares and P. Brax, _Detecting coupled domain walls in laboratory experiments_ , _Phys. Rev. Lett._ 122 (2019) 091102.
  * [25] Y.-C. Cai, N. Padilla and B. Li, _Testing gravity using cosmic voids_ , _MNRAS_ 451 (2015) 1036 [1410.1510].
  * [26] E.L.D. Perico, R. Voivodic, M. Lima and D.F. Mota, _Cosmic voids in modified gravity scenarios_ , _arXiv e-prints_ (2019) arXiv:1905.12450 [1905.12450].
  * [27] N.D. Padilla, _Testing Gravity with Cosmic Voids_ , _Boletin de la Asociacion Argentina de Astronomia La Plata Argentina_ 57 (2015) 6.
  * [28] C.T. Davies, M. Cautun and B. Li, _Cosmological test of gravity using weak lensing voids_ , _MNRAS_ 490 (2019) 4907 [1907.06657].
  * [29] J. Clampitt, Y.-C. Cai and B. Li, _Voids in modified gravity: excursion set predictions_ , _MNRAS_ 431 (2013) 749 [https://academic.oup.com/mnras/article-pdf/431/1/749/18243874/stt219.pdf].
  * [30] T. Chantavat, U. Sawangwit and B.D. Wandelt, _Void Profile from Planck Lensing Potential Map_ , _ApJ_ 836 (2017) 156 [1702.01009].
  * [31] W. Hu and I. Sawicki, _Models of f(R) cosmic acceleration that evade solar system tests_ , _Phys. Rev. D_ 76 (2007) 064004 [0705.1158].
  * [32] M. Turner, R. Clough, H. Martin and L. Topp, _Stiffness and deflection analysis of complex structures_ , _Journal of Aeronautical Sciences_ 23 (1956) 805.
  * [33] W.K. Liu, S. Li and H. Park, _Eighty Years of the Finite Element Method: Birth, Evolution, and Future_ , _arXiv e-prints_ (2021) arXiv:2107.04960 [2107.04960].
  * [34] A. Hrennikoff, _Solution of Problems of Elasticity by the Framework Method_ , _Journal of Applied Mechanics_ 8 (2021) A169.
  * [35] J. Braden, C. Burrage, B. Elder and D. Saadeh, _$\varphi$ enics: Vainshtein screening with the finite element method_, _J. Cosmology Astropart. Phys_ 2021 (2021) 010 [2011.07037].
  * [36] C. Burrage, B. Elder, P. Millington, D. Saadeh and B. Thrussell, _Fifth-force screening around extremely compact sources_ , _J. Cosmology Astropart. Phys_ 2021 (2021) 052 [2104.14564].
  * [37] C. Briddon, C. Burrage, A. Moss and A. Tamosiunas, _SELCIE: A tool for investigating the chameleon field of arbitrary sources_ , _arXiv e-prints_ (2021) arXiv:2110.11917 [2110.11917].
  * [38] A. Logg and G.N. Wells, _DOLFIN: Automated Finite Element Computing_ , _arXiv e-prints_ (2011) arXiv:1103.6248 [1103.6248].
  * [39] A. Logg, K.-A. Mardal, G.N. Wells et al., _Automated Solution of Differential Equations by the Finite Element Method_ , Springer (2012), 10.1007/978-3-642-23099-8.
  * [40] M.S. Alnaes, A. Logg, K.B. Oelgaard, M.E. Rognes and G.N. Wells, _Unified Form Language: A domain-specific language for weak formulations of partial differential equations_ , _arXiv e-prints_ (2012) arXiv:1211.4047 [1211.4047].
  * [41] M.W. Scroggs, J.S. Dokken, C.N. Richardson and G.N. Wells, _Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes_ , _arXiv e-prints_ (2021) arXiv:2102.11901 [2102.11901].
  * [42] M. Plionis and S. Basilakos, _The size and shape of local voids_ , _MNRAS_ 330 (2002) 399 [astro-ph/0106491].
  * [43] F. Hoyle and M.S. Vogeley, _Voids in the Point Source Catalogue Survey and the Updated Zwicky Catalog_ , _ApJ_ 566 (2002) 641 [astro-ph/0109357].
  * [44] E. Russell and J.R. Pycke, _Log-Normal Distribution of Cosmic Voids in Simulations and Mocks_ , _ApJ_ 835 (2017) 69 [1612.03180].
  * [45] R.P. Kirshner, J. Oemler, Augustus, P.L. Schechter and S.A. Shectman, _A Survey of the Bootes Void_ , _ApJ_ 314 (1987) 493.
  * [46] G.A. Wegner, J.J. Salzer, J.M. Taylor and A.S. Hirschauer, _Metal Abundances and Star Formation Rates of Emission-line Galaxies in and around the Boötes Void_ , _ApJ_ 883 (2019) 29 [1908.07539].
  * [47] I.D. Karachentsev, V.E. Karachentseva, W.K. Huchtmeier and D.I. Makarov, _A catalog of neighboring galaxies_ , _AJ_ 127 (2004) 2031.
  * [48] R. van de Weygaert and E. Platen, _Cosmic Voids: Structure, Dynamics and Galaxies_ , in _International Journal of Modern Physics Conference Series_ , vol. 1 of _International Journal of Modern Physics Conference Series_ , pp. 41–66, Jan., 2011, DOI [0912.2997].
  * [49] N. Hamaus, P.M. Sutter and B.D. Wandelt, _Universal Density Profile for Cosmic Voids_ , _Phys. Rev. Lett._ 112 (2014) 251302 [1403.5499].
  * [50] S. Nadathur, S. Hotchkiss, J.M. Diego, I.T. Iliev, S. Gottlöber, W.A. Watson et al., _Universal void density profiles from simulation and SDSS_ , in _The Zeldovich Universe: Genesis and Growth of the Cosmic Web_ , R. van de Weygaert, S. Shandarin, E. Saar and J. Einasto, eds., vol. 308, pp. 542–545, Oct., 2016, DOI [1412.8372].
  * [51] S. Nadathur, S. Hotchkiss, J.M. Diego, I.T. Iliev, S. Gottlöber, W.A. Watson et al., _Self-similarity and universality of void density profiles in simulation and SDSS data_ , _MNRAS_ 449 (2015) 3997 [https://academic.oup.com/mnras/article-pdf/449/4/3997/18505550/stv513.pdf].
  * [52] S. Nadathur and S. Hotchkiss, _A robust public catalogue of voids and superclusters in the SDSS Data Release 7 galaxy surveys_ , _MNRAS_ 440 (2014) 1248 [1310.2791].
  * [53] W.A. Watson, J.M. Diego, S. Gottlöber, I.T. Iliev, A. Knebe, E. Martínez-González et al., _The Jubilee ISW project – I. Simulated ISW and weak lensing maps and initial power spectra results_ , _MNRAS_ 438 (2013) 412 [https://academic.oup.com/mnras/article-pdf/438/1/412/18462721/stt2208.pdf].
  * [54] M.C. Neyrinck, _zobov: a parameter-free void-finding algorithm_ , _MNRAS_ 386 (2008) 2101 [https://academic.oup.com/mnras/article-pdf/386/4/2101/4276375/mnras0386-2101.pdf].
  * [55] E. Ricciardelli, V. Quilis and J. Varela, _On the universality of void density profiles_ , _MNRAS_ 440 (2014) 601 [1402.2976].
  * [56] E. Massara and R.K. Sheth, _Density and velocity profiles around cosmic voids_ , _arXiv e-prints_ (2018) arXiv:1811.03132 [1811.03132].
  * [57] J. Velasquez and L. Castaneda, _Equivalence between Scalar-Tensor theories and $f(R)$-gravity: From the action to Cosmological Perturbations_, _arXiv e-prints_ (2018) arXiv:1808.05615 [1808.05615].
  * [58] L. Pizzuti, B. Sartoris, L. Amendola, S. Borgani, A. Biviano, K. Umetsu et al., _Clash-vlt: constraints on f(r) gravity models with galaxy clusters using lensing and kinematic analyses_ , _Journal of Cosmology and Astroparticle Physics_ 2017 (2017) 023.
  * [59] S. Shao, B. Li, M. Cautun, H. Wang and J. Wang, _Screening maps of the local Universe I - Methodology_ , _MNRAS_ 489 (2019) 4912 [1907.02081].
  * [60] C. Wilson and R. Bean, _Testing f (R ) gravity with scale dependent cosmic void velocity profiles_ , _Phys. Rev. D_ 104 (2021) 023512 [2012.05925].