
Type Ia supernovae (SNIe) are generally accepted to act as standardisable
candles, and their use in cosmology led to the first confirmation of the as yet
unexplained accelerated cosmic expansion. Many of the theoretical models to
explain the cosmic acceleration assume modifications to Einsteinian General
Relativity which accelerate the expansion, but the question of whether such
modifications also affect the ability of SNIe to be standardisable candles has
rarely been addressed. This paper is an attempt to answer this question. For
this we adopt a semianalytical model to calculate SNIe light curves in
nonstandard gravity. We use this model to show that the average rescaled
intrinsic peak luminosity  a quantity that is assumed to be constant with
redshift in standard analyses of Type Ia supernova (SNIa) cosmology data 
depends on the strength of gravity in the supernova's local environment because
the latter determines the Chandrasekhar mass  the mass of the SNIa's white
dwarf progenitor right before the explosion. This means that SNIe are no longer
standardisable candles in scenarios where the strength of gravity evolves over
time, and therefore the cosmology implied by the existing SNIa data will be
different when analysed in the context of such models. As an example, we show
that the observational SNIa cosmology data can be fitted with both a model
where $(\Omega_{\rm M}, \Omega_{\Lambda})=(0.62, 0.38)$ and Newton's constant
$G$ varies as $G(z)=G_0(1+z)^{1/4}$ and the standard model where $(\Omega_{\rm
M}, \Omega_{\Lambda})=(0.3, 0.7)$ and $G$ is constant, when the Universe is
assumed to be flat.

The intrinsic peak luminosity of Type Ia supernovae (SNe Ia) depends on the
value of Newton's gravitational constant $G$, through the Chandrasekhar mass
$M_{\rm Ch}\propto G^{3/2}$. If the luminosity distance $d_{\rm L}$ can be
independently determined, the SNe Ia can be treated as a tracker to constrain
the possible time variation of $G$ in different redshift ranges. The
gravitationalwave (GW) standard sirens, caused by the coalescence of binary
neutron stars, provide a modelindependent way to measure the distance of GW
events, which can be used to determine the luminosity distances of SNe Ia by
interpolation, provided the GW and SN Ia samples have similar redshift ranges.
We demonstrate that combining the GW observations of thirdgeneration detectors
with SN Ia data provides a powerful and modelindependent way to measure $G$ in
a wide redshift range, which can constrain the ratio $G/G_0$, where $G$ and
$G_0$ are respectively the values in the redshift ranges $z>0.1$ and $z<0.1$,
at the level of $1.5\%$.

We propose a new filter, a smooth$k$ space filter, to use in the
PressSchechter approach to model the dark matter halo mass function which
overcomes shortcomings of other filters. We test this against the mass function
measured in Nbody simulations. We find that the commonly used sharp$k$ filter
fails to reproduce the behaviour of the halo mass function at low masses
measured from simulations of models with a sharp truncation in the linear power
spectrum. We show that the predictions with our new filter agree with the
simulation results over a wider range of halo masses for both damped and
undamped power spectra than is the case with the sharp$k$ and realspace
tophat filters.

Cosmic voids are an important probe of largescale structure that allows us
to constrain cosmological parameters and test cosmological models. We present a
new paradigm for void studies: void detection in weak lensing convergence maps.
This approach identifies objects that relate directly to our theoretical
understanding of voids as underdensities in the total matter field and presents
several advantages compared to the customary method of finding voids in the
galaxy distribution. We exemplify this approach by identifying voids using the
weak lensing peaks as tracers of the largescale structure. We find
selfsimilarity in the void abundance across a range of peak signaltonoise
selection thresholds. The voids obtained via this approach give a tangential
shear signal up to $\sim50$ times larger than voids identified in the galaxy
distribution.

We present a systematic comparison of several existing and new void finding
algorithms, focusing on their potential power to test a particular class of
modified gravity models  chameleon $f(R)$ gravity. These models deviate from
standard General Relativity (GR) more strongly in lowdensity regions and thus
voids are a promising venue to test them. We use Halo Occupation Distribution
(HOD) prescriptions to populate haloes with galaxies, and tune the HOD
parameters such that the galaxy twopoint correlation functions are the same in
both f(R) and GR models. We identify both 3D voids as well as 2D underdensities
in the planeofthesky to find the same void abundance and void galaxy number
density profiles across all models, which suggests that they do not contain
much information beyond galaxy clustering. However, the underlying void dark
matter density profiles are significantly different, with f(R) voids being more
underdense than GR ones, which leads to f(R) voids having a larger tangential
shear signal than their GR analogues. We investigate the potential of each void
finder to test f(R) models with nearfuture lensing surveys such as EUCLID and
LSST. The 2D voids have the largest power to probe f(R) gravity, with a LSST
analysis of tunnel (which is a new type of 2D underdensity introduced here)
lensing distinguishing at 80 and 11$\sigma$ (statistical error) f(R) models
with $f_{R0}=10^{5}$ and $10^{6}$ from GR.

We propose a new framework for testing gravity using cluster observations,
which aims to provide an unbiased constraint on modified gravity models from
Sunyaev Zel'dovich (SZ) and Xray cluster counts and the cluster gas fraction,
among other possible observables. Focusing on a popular $ f(R) $ model of
gravity, we propose a novel procedure to recalibrate mass scaling relations
from $ \Lambda $CDM to $ f(R) $ gravity for SZ and Xray cluster observables.
We find that the complicated modified gravity effects can be simply modelled as
a dependence on a combination of the background scalar field and redshift, $
f_R(z)/(1+z) $, regardless of the $ f(R) $ model parameter. By employing a
large suite of Nbody simulations, we demonstrate that a theoretically derived
tanh fitting formula is in excellent agreement with the dynamical mass
enhancement of dark matter haloes for a large range of background field
parameters and redshifts. Our framework is sufficiently flexible to allow for
tests of other models and inclusion of further observables. The oneparameter
description of the dynamical mass enhancement can have important implications
on the theoretical modelling of observables and on practical tests of gravity.

In theories of modified gravity with the chameleon screening mechanism, the
strength of the fifth force depends on environment. This induces an environment
dependence of structure formation, which differs from $\Lambda$CDM. We show
that these differences can be captured by the marked correlation function. With
the galaxy correlation functions and number densities calibrated to match
between $f(R)$ and $\Lambda$CDM models in simulations, we show that the marked
correlation functions from using either the local density or halo mass as the
marks encode extra information, which can be used to test these theories. We
discuss possible applications of these statistics in observations.

We analyse the twopoint and marked correlation functions of haloes and
galaxies in three variants of the chameleon $f(R)$ gravity model using Nbody
simulations, and compare to a fiducial $\Lambda$CDM model based on general
relativity (GR). Using a halo occupation distribution prescription (HOD) we
populate dark matter haloes with galaxies, where the HOD parameters have been
tuned such that the galaxy number densities and the realspace galaxy twopoint
correlation functions in the modified gravity models match those in GR to
within $1\sim3\%$. We test the idea that since the behaviour of gravity is
dependent on environment, marked correlation functions may display a measurable
difference between the models. For this we test marks based on the density
field and the Newtonian gravitational potential. We find that the galaxy marked
correlation function shows significant differences measured in different models
on scales smaller than $r\lesssim 20~h^{1}$ Mpc. Guided by simulations to
identify a suitable mark, this approach could be used as a new probe of the
accelerated expansion of the Universe.

A theoretically interesting and practically important question in cosmology
is the reconstruction of the initial density distribution provided a latetime
density field. This is a longstanding question with a revived interest
recently, especially in the context of optimally extracting the baryonic
acoustic oscillation (BAO) signals from observed galaxy distributions. We
present a new efficient method to carry out this reconstruction, which is based
on numerical solutions to the nonlinear partial differential equation that
governs the mapping between the initial Lagrangian and final Eulerian
coordinates of particles in evolved density fields. This is motivated by
numerical simulations of the quartic Galileon gravity model, which has similar
equations that can be solved effectively by multigrid GaussSeidel relaxation.
The method is based on mass conservation, and does not assume any specific
cosmological model. Our test shows that it has a performance comparable to that
of stateoftheart algorithms which were very recently put forward in the
literature, with the reconstructed density field over $\sim80\%$ ($50\%$)
correlated with the initial condition at $k\lesssim0.6h/{\rm Mpc}$ ($1.0h/{\rm
Mpc}$). With an example, we demonstrate that this method can significantly
improve the accuracy of BAO reconstruction.

We present an analysis of galaxygalaxy weak gravitational lensing (GGL) in
chameleon $f(R)$ gravity  a leading candidate of nonstandard gravity models.
For the analysis we have created mock galaxy catalogues based on dark matter
haloes from two sets of numerical simulations, using a halo occupation
distribution (HOD) prescription which allows a redshift dependence of galaxy
number density. To make a fairer comparison between the $f(R)$ and $\Lambda$CDM
models, their HOD parameters are tuned so that the galaxy twopoint correlation
functions in real space (and therefore the projected twopoint correlation
functions) match. While the $f(R)$ model predicts an enhancement of the
convergence power spectrum by up to $\sim30\%$ compared to the standard
$\Lambda$CDM model with the same parameters, the maximum enhancement of GGL is
only half as large and less than 5\% on separations above $\sim1$$2h^{1}$Mpc,
because the latter is a cross correlation of shear (or matter, which is more
strongly affected by modified gravity) and galaxy (which is weakly affected
given the good match between galaxy auto correlations in the two models)
fields. We also study the possibility of reconstructing the matter power
spectrum by combination of GGL and galaxy clustering in $f(R)$ gravity. We find
that the galaxymatter cross correlation coefficient remains at unity down to
$\sim2$$3h^{1}$Mpc at relevant redshifts even in $f(R)$ gravity, indicating
joint analysis of GGL and galaxy clustering can be a powerful probe of matter
density fluctuations in chameleon gravity. The scale dependence of the model
differences in their predictions of GGL can potentially allow to break the
degeneracy between $f(R)$ gravity and other cosmological parameters such as
$\Omega_m$ and $\sigma_8$.

Scalar tensor theories can be expressed in different frames, such as the
commonlyused Einstein and Jordan frames, and it is generally accepted that
cosmological observables are the same in these frames. We revisit this by
making a detailed sidebyside comparison of the quantities and equations in
two conformally related frames, from the actions and fully covariant field
equations to the linearised equations in both real and Fourier spaces. This
confirms that the field and conservation equations are equivalent in the two
frames, in the sense that we can always reexpress equations in one frame using
relevant transformations of variables to derive the corresponding equations in
the other. We show, with both analytical derivation and a numerical example,
that the lineofsight integration to calculate CMB temperature anisotropies
can be done using either Einstein frame or Jordan frame quantities, and the
results are identical, provided the correct redshift is used in the Einstein
frame ($1+z\neq1/a$).

We investigate the role of thermal velocities in Nbody simulations of
structure formation in warm dark matter models. Starting from the commonly used
approach of adding thermal velocities, randomly selected from a FermiDirac
distribution, to the gravitationallyinduced (peculiar) velocities of the
simulation particles, we compare the matter and velocity power spectra measured
from CDM and WDM simulations with and without thermal velocities. This
prescription for adding thermal velocities results in deviations in the
velocity field in the initial conditions away from the linear theory
predictions, which affects the evolution of structure at later times. We show
that this is entirely due to numerical noise. For a warm candidate with mass
$3.3$ keV, the matter and velocity power spectra measured from simulations with
thermal velocities starting at $z=199$ deviate from the linear prediction at $k
\gtrsim10$ $h/$Mpc, with an enhancement of the matter power spectrum $\sim
\mathcal{O}(10)$ and of the velocity power spectrum $\sim \mathcal{O}(10^2)$ at
wavenumbers $k \sim 64$ $h/$Mpc with respect to the case without thermal
velocities. At late times, these effects tend to be less pronounced. Indeed, at
$z=0$ the deviations do not exceed $6\%$ (in the velocity spectrum) and $1\%$
(in the matter spectrum) for scales $10 <k< 64$ $h/$Mpc. Increasing the
resolution of the Nbody simulations shifts these deviations to higher
wavenumbers. The noise introduces more spurious structures in WDM simulations
with thermal velocities and modifies the radial density profiles of dark matter
haloes. We find that spurious haloes start to appear in simulations which
include thermal velocities at a mass that is $\sim$3 times larger than in
simulations without thermal velocities.

We refine the mass and environment dependent spherical collapse model of
chameleon $f(R)$ gravity by calibrating a phenomenological correction inspired
by the parameterized postFriedmann framework against highresolution $N$body
simulations. We employ our method to predict the corresponding modified halo
mass function, and provide fitting formulas to calculate the fractional
enhancement of the $f(R)$ halo abundance with respect to that of General
Relativity (GR) within a precision of $\lesssim 5\%$ from the results obtained
in the simulations. Similar accuracy can be achieved for the full $f(R)$ mass
function on the condition that the modeling of the reference GR abundance of
halos is accurate at the percent level. We use our fits to forecast constraints
on the additional scalar degree of freedom of the theory, finding that upper
bounds competitive with current Solar System tests are within reach of cluster
number count analyses from ongoing and upcoming surveys at much larger scales.
Importantly, the flexibility of our method allows also for this to be applied
to other scalartensor theories characterized by a mass and environment
dependent spherical collapse.

We describe and demonstrate the potential of a new and very efficient method
for simulating certain classes of modified gravity theories, such as the widely
studied $f(R)$ gravity models. High resolution simulations for such models are
currently very slow due to the highly nonlinear partial differential equation
that needs to be solved exactly to predict the modified gravitational force.
This nonlinearity is partly inherent, but is also exacerbated by the specific
numerical algorithm used, which employs a variable redefinition to prevent
numerical instabilities. The standard NewtonGaussSeidel iterative method used
to tackle this problem has a poor convergence rate. Our new method not only
avoids this, but also allows the discretised equation to be written in a form
that is analytically solvable. We show that this new method greatly improves
the performance and efficiency of $f(R)$ simulations. For example, a test
simulation with $512^3$ particles in a box of size $512 \, \mathrm{Mpc}/h$ is
now 5 times faster than before, while a Millenniumresolution simulation for
$f(R)$ gravity is estimated to be more than 20 times faster than with the old
method. Our new implementation will be particularly useful for running very
high resolution, largesized simulations which, to date, are only possible for
the standard model, and also makes it feasible to run large numbers of lower
resolution simulations for covariance analyses. We hope that the method will
bring us to a new era for precision cosmological tests of gravity.

We investigate the information content of various cosmic shear statistics on
the theory of gravity. Focusing on the HuSawickitype $f(R)$ model, we perform
a set of raytracing simulations and measure the convergence bispectrum, peak
counts and Minkowski functionals. We first show that while the convergence
power spectrum does have sensitivity to the current value of extra scalar
degree of freedom $f_{\rm R0}$, it is largely compensated by a change in the
present density amplitude parameter $\sigma_{8}$ and the matter density
parameter $\Omega_{\rm m0}$. With accurate covariance matrices obtained from
1000 lensing simulations, we then examine the constraining power of the three
additional statistics. We find that these probes are indeed helpful to break
the parameter degeneracy, which can not be resolved from the power spectrum
alone. We show that especially the peak counts and Minkowski functionals have
the potential to rigorously (marginally) detect the signature of modified
gravity with the parameter $f_{\rm R0}$ as small as $10^{5}$ ($10^{6}$) if
we can properly model them on small ($\sim 1\, \mathrm{arcmin}$) scale in a
future survey with a sky coverage of 1,500 squared degrees. We also show that
the signal level is similar among the additional three statistics and all of
them provide complementary information to the power spectrum. These findings
indicate the importance of combining multiple probes beyond the standard power
spectrum analysis to detect possible modifications to General Relativity.

DESI (Dark Energy Spectroscopic Instrument) is a Stage IV groundbased dark
energy experiment that will study baryon acoustic oscillations (BAO) and the
growth of structure through redshiftspace distortions with a widearea galaxy
and quasar redshift survey. To trace the underlying dark matter distribution,
spectroscopic targets will be selected in four classes from imaging data. We
will measure luminous red galaxies up to $z=1.0$. To probe the Universe out to
even higher redshift, DESI will target bright [O II] emission line galaxies up
to $z=1.7$. Quasars will be targeted both as direct tracers of the underlying
dark matter distribution and, at higher redshifts ($ 2.1 < z < 3.5$), for the
Ly$\alpha$ forest absorption features in their spectra, which will be used to
trace the distribution of neutral hydrogen. When moonlight prevents efficient
observations of the faint targets of the baseline survey, DESI will conduct a
magnitudelimited Bright Galaxy Survey comprising approximately 10 million
galaxies with a median $z\approx 0.2$. In total, more than 30 million galaxy
and quasar redshifts will be obtained to measure the BAO feature and determine
the matter power spectrum, including redshift space distortions.

DESI (Dark Energy Spectropic Instrument) is a Stage IV groundbased dark
energy experiment that will study baryon acoustic oscillations and the growth
of structure through redshiftspace distortions with a widearea galaxy and
quasar redshift survey. The DESI instrument is a roboticallyactuated,
fiberfed spectrograph capable of taking up to 5,000 simultaneous spectra over
a wavelength range from 360 nm to 980 nm. The fibers feed ten threearm
spectrographs with resolution $R= \lambda/\Delta\lambda$ between 2000 and 5500,
depending on wavelength. The DESI instrument will be used to conduct a
fiveyear survey designed to cover 14,000 deg$^2$. This powerful instrument
will be installed at prime focus on the 4m Mayall telescope in Kitt Peak,
Arizona, along with a new optical corrector, which will provide a threedegree
diameter field of view. The DESI collaboration will also deliver a
spectroscopic pipeline and data management system to reduce and archive all
data for eventual public use.

Abell 2146 ($z$ = 0.232) consists of two galaxy clusters undergoing a major
merger. The system was discovered in previous work, where two large shock
fronts were detected using the $\textit{Chandra Xray Observatory}$, consistent
with a merger close to the plane of the sky, caught soon after first core
passage. A weak gravitational lensing analysis of the total gravitating mass in
the system, using the distorted shapes of distant galaxies seen with ACSWFC on
$\textit{Hubble Space Telescope}$, is presented. The highest peak in the
reconstruction of the projected mass is centred on the Brightest Cluster Galaxy
(BCG) in Abell 2146A. The mass associated with Abell 2146B is more extended.
Bootstrapped noise mass reconstructions show the mass peak in Abell 2146A to
be consistently centred on the BCG. Previous work showed that BCGA appears to
lag behind an Xray cool core; although the peak of the mass reconstruction is
centred on the BCG, it is also consistent with the Xray peak given the
resolution of the weak lensing mass map. The bestfit mass model with two
components centred on the BCGs yields $M_{200}$ =
1.1$^{+0.3}_{0.4}$$\times$10$^{15}$M$_{\odot}$ and
3$^{+1}_{2}$$\times$10$^{14}$M$_{\odot}$ for Abell 2146A and Abell 2146B
respectively, assuming a mass concentration parameter of $c=3.5$ for each
cluster. From the weak lensing analysis, Abell 2146A is the primary halo
component, and the origin of the apparent discrepancy with the Xray analysis
where Abell 2146B is the primary halo is being assessed using simulations of
the merger.

In this Letter, we report the observational constraints on the HuSawicki
$f(R)$ theory derived from weak lensing peak abundances, which are closely
related to the mass function of massive halos. In comparison with studies using
optical or xray clusters of galaxies, weak lensing peak analyses have the
advantages of not relying on massbaryonic observable calibrations. With
observations from the CanadaFranceHawaiiTelescope Lensing Survey, our peak
analyses give rise to a tight constraint on the model parameter $f_{R0}$ for
$n=1$. The $95\%$ CL limit is $\log_{10}f_{R0} < 4.82$ given WMAP9 priors on
$(\Omega_{\rm m}, A_{\rm s})$. With Planck15 priors, the corresponding result
is $\log_{10}f_{R0} < 5.16$.

Modified theories of gravity provide us with a unique opportunity to generate
innovative tests of gravity. In Chameleon f(R) gravity, the gravitational
potential differs from the weakfield limit of general relativity (GR) in a
mass dependent way. We develop a probe of gravity which compares high mass
clusters, where Chameleon effects are weak, to low mass clusters, where the
effects can be strong. We utilize the escape velocity edges in the
radius/velocity phase space to infer the gravitational potential profiles on
scales of 0.31 virial radii. We show that the escape edges of low mass
clusters are enhanced compared to GR, where the magnitude of the difference
depends on the background field value fR0. We validate our probe using Nbody
simulations and simulated light cone galaxy data. For a DESI (Dark Energy
Spectroscopic Instrument) Bright Galaxy Sample, including observational
systematics, projection effects, and cosmic variance, our test can
differentiate between GR and Chameleon f(R) gravity models, fR0 = 4e6 (2e6)
at > 5{\sigma} (> 2{\sigma}), more than an order of magnitude better than
current clusterscale constraints.

We introduce and demonstrate the power of a method to speed up current
iterative techniques for Nbody modified gravity simulations. Our method is
based on the observation that the accuracy of the final result is not
compromised if the calculation of the fifth force becomes less accurate, but
substantially faster, in highdensity regions where it is weak due to
screening. We focus on the nDGP model which employs Vainshtein screening, and
test our method by running AMR simulations in which the solutions on the finer
levels of the mesh (high density) are not obtained iteratively, but instead
interpolated from coarser levels. We show that the impact this has on the
matter power spectrum is below $1\%$ for $k < 5h/{\rm Mpc}$ at $z = 0$, and
even smaller at higher redshift. The impact on halo properties is also small
($\lesssim 3\%$ for abundance, profiles, mass; and $\lesssim 0.05\%$ for
positions and velocities). The method can boost the performance of modified
gravity simulations by more than a factor of 10, which allows them to be pushed
to resolution levels that were previously hard to achieve.

Many modified gravity theories are under consideration in cosmology as the
source of the accelerated expansion of the universe and linear perturbation
theory, valid on the largest scales, has been examined in many of these models.
However, smaller nonlinear scales offer a richer phenomenology with which to
constrain modified gravity theories. Here, we consider the HuSawicki form of
$f(R)$ gravity and apply the postFriedmann approach to derive the leading
order equations for nonlinear scales, i.e. the equations valid in the
Newtonianlike regime. We reproduce the standard equations for the scalar
field, gravitational slip and the modified Poisson equation in a coherent
framework. In addition, we derive the equation for the leading order correction
to the Newtonian regime, the vector potential. We measure this vector potential
from $f(R)$ Nbody simulations at redshift zero and one, for two values of the
$f_{R_0}$ parameter. We find that the vector potential at redshift zero in
$f(R)$ gravity can be close to 50\% larger than in GR on small scales for
$f_{R_0}=1.289\times10^{5}$, although this is less for larger scales,
earlier times and smaller values of the $f_{R_0}$ parameter. Similarly to in
GR, the small amplitude of this vector potential suggests that the Newtonian
approximation is highly accurate for $f(R)$ gravity, and also that the
nonlinear cosmological behaviour of $f(R)$ gravity can be completely described
by just the scalar potentials and the $f(R)$ field.

Selfconsistent ${\it N}$body simulations of modified gravity models are a
key ingredient to obtain rigorous constraints on deviations from General
Relativity using largescale structure observations. This paper provides the
first detailed comparison of the results of different ${\it N}$body codes for
the $f(R)$, DGP, and Symmetron models, starting from the same initial
conditions. We find that the fractional deviation of the matter power spectrum
from $\Lambda$CDM agrees to better than $1\%$ up to $k \sim 510~h/{\rm Mpc}$
between the different codes. These codes are thus able to meet the stringent
accuracy requirements of upcoming observational surveys. All codes are also in
good agreement in their results for the velocity divergence power spectrum,
halo abundances and halo profiles. We also test the quasistatic limit, which
is employed in most modified gravity ${\it N}$body codes, for the Symmetron
model for which the most significant nonstatic effects among the models
considered are expected. We conclude that this limit is a very good
approximation for all of the observables considered here.

We explore the Minkowski functionals of weak lensing convergence map to
distinguish between $f(R)$ gravity and the general relativity (GR). The mock
weak lensing convergence maps are constructed with a set of highresolution
simulations assuming different gravity models. It is shown that the lensing MFs
of $f(R)$ gravity can be considerably different from that of GR because of the
environmentally dependent enhancement of structure formation. We also
investigate the effect of lensing noise on our results, and find that it is
likely to distinguish F5, F6 and GR gravity models with a galaxy survey of
$\sim3000$ degree$^2$ and with a background source number density of
$n_g=30~{\rm arcmin}^{2}$, comparable to an upcoming survey dark energy survey
(DES). We also find that the $f(R)$ signal can be partially degenerate with the
effect of changing cosmology, but combined use of other observations, such as
the cosmic microwave background (CMB) data, can help break this degeneracy.

We study lensing by voids in Cubic Galileon and Nonlocal gravity cosmologies,
which are examples of theories of gravity that modify the lensing potential. We
find voids in the dark matter and halo density fields of Nbody simulations and
compute their lensing signal analytically from the void density profiles, which
we show are well fit by a simple analytical formula. In the Cubic Galileon
model, the modifications to gravity inside voids are not screened and they
approximately double the size of the lensing effects compared to GR. The
difference is largely determined by the direct effects of the fifth force on
lensing and less so by the modified density profiles. For this model, we also
discuss the subtle impact on the force and lensing calculations caused by the
screening effects of haloes that exist in and around voids. In the Nonlocal
model, the impact of the modified density profiles and the direct modifications
to lensing are comparable, but they boost the lensing signal by only $\approx
10\%$, compared with that of GR. Overall, our results suggest that lensing by
voids is a promising tool to test models of gravity that modify lensing.