We investigate the accuracy of weak lensing simulations by comparing the results of five independently developed lensing simulation codes run on the same input $N$-body simulation. Our comparison focuses on the lensing convergence maps produced by the codes, and in particular on the corresponding PDFs, power spectra and peak counts. We find that the convergence power spectra of the lensing codes agree to $\lesssim 2\%$ out to scales $\ell \approx 4000$. For lensing peak counts, the agreement is better than $5\%$ for peaks with signal-to-noise $\lesssim 6$. We also discuss the systematic errors due to the Born approximation, line-of-sight discretization, particle noise and smoothing. The lensing codes tested deal in markedly different ways with these effects, but they nonetheless display a satisfactory level of agreement. Our results thus suggest that systematic errors due to the operation of existing lensing codes should be small. Moreover their impact on the convergence power spectra for a lensing simulation can be predicted given its numerical details, which may then serve as a validation test.

We study the impact that baryon-CDM relative density perturbations $\delta_{bc}$ have on galaxy formation using cosmological simulations with the IllustrisTNG model. These isocurvature (non-adiabatic) perturbations can be induced primordially, if multiple fields are present during inflation, and are generated before baryon-photon decoupling when baryons did not comove with CDM. The presence of long-wavelength $\delta_{bc}$ perturbations in our simulations is mimicked by modifying the ratios of the cosmic densities of baryons $\Omega_b$ and CDM $\Omega_c$, at fixed total matter density $\Omega_m$. We measure the corresponding galaxy bias parameter $b_{\delta}^{bc}$ as the response of galaxy abundances to $\delta_{bc}$. When selecting by total host halo mass, $b_{\delta}^{bc}$ is negative and it decreases with mass and redshift. Stellar-mass selected simulated galaxies show a weaker or even the opposite trend because of the competing effects of $\delta_{bc}$ on the halo mass function and stellar-to-halo-mass relations. We show that simple modeling of the latter two effects describes $b_{\delta}^{bc}$ for stellar-mass-selected objects well. We find $b_{\delta}^{bc} =0.6$ for $M_* = 10^{11}\ M_{\odot}/h$ and $z=0.5$, which is representative of BOSS DR12 galaxies. For $\delta_{bc}$ modes generated by baryon-photon interactions, we estimate the impact on the DR12 power spectrum to be below $1\%$, and shifts on inferred distance and growth rate parameters should not exceed $0.1\%$.

We measure power spectrum response functions in the presence of baryonic physical processes using separate universe simulations with the IllustrisTNG galaxy formation model. The response functions describe how the small-scale power spectrum reacts to long-wavelength perturbations and they can be efficiently measured with the separate universe technique by absorbing the effects of the long modes into a modified cosmology. Specifically, we focus on the total first-order matter power spectrum response to an isotropic density fluctuation $R_1(k,z)$, which is fully determined by the logarithmic derivative of the nonlinear matter power spectrum ${\rm dln}P_m(k,z)/{\rm dln}k$ and the growth-only response function $G_1(k,z)$. We find that $G_1(k,z)$ is not affected by the baryonic physical processes in the simulations at redshifts $z < 3$ and on all scales probed ($k \lesssim 15h/{\rm Mpc}$, i.e. length scales $\gtrsim 0.4 {\rm Mpc}/h$). In practice, this implies that the power spectrum fully specifies the baryonic dependence of its response function. Assuming an idealized lensing survey setup, we evaluate numerically the baryonic impact on the squeezed-lensing bispectrum and the lensing super-sample power spectrum covariance, which are given in terms of responses. Our results show that these higher-order lensing statistics can display varying levels of sensitivity to baryonic effects compared to the power spectrum, with the squeezed-bispectrum being the least sensitive. We also show that ignoring baryonic effects on lensing covariances slightly overestimates the error budget (and is therefore conservative from the point of view of parameter error bars) and likely has negligible impact on parameter biases in inference analyses.

We study the impact that large-scale perturbations of (i) the matter density and (ii) the primordial gravitational potential with local primordial non-Gaussianity (PNG) have on galaxy formation using the IllustrisTNG model. We focus on the linear galaxy bias $b_1$ and the coefficient $b_\phi$ of the scale-dependent bias induced by PNG, which describe the response of galaxy number counts to these two types of perturbations, respectively. We perform our study using separate universe simulations, in which the effect of the perturbations is mimicked by changes to the cosmological parameters: modified cosmic matter density for $b_1$ and modified amplitude $\mathcal{A}_s$ of the primordial scalar power spectrum for $b_\phi$. We find that the widely used universality relation $b_\phi = 2\delta_c(b_1 - 1)$ is a poor description of the bias of haloes and galaxies selected by stellar mass $M_*$, which is instead described better by $b_\phi(M_*) = 2\delta_c(b_1(M_*) - p)$ with $p \in [0.4, 0.7]$. This is explained by the different impact that matter overdensities and local PNG have on the median stellar-to-halo-mass relation. A simple model of this impact allows us to describe the stellar mass dependence of $b_1$ and $b_\phi$ fairly well. Our results also show a nontrivial relation between $b_1$ and $b_\phi$ for galaxies selected by color and black hole mass accretion rate. Our results provide refined priors on $b_\phi$ for local PNG constraints and forecasts using galaxy clustering. Given that the widely used universality relation underpredicts $b_\phi(M_*)$, existing analyses may underestimate the true constraining power on local PNG.

We investigate the potential of the galaxy power spectrum to constrain compensated isocurvature perturbations (CIPs), primordial fluctuations in the baryon density that are compensated by fluctuations in CDM density to ensure an unperturbed total matter density. We show that CIPs contribute to the galaxy overdensity at linear order, and if they are close to scale-invariant, their effects are nearly perfectly degenerate with the local PNG parameter $f_{\rm nl}$ if they correlate with the adiabatic perturbations. This degeneracy can however be broken by analyzing multiple galaxy samples with different bias parameters, or by taking CMB priors on $f_{\rm nl}$ into account. Parametrizing the amplitude of the CIP power spectrum as $P_{\sigma\sigma} = A_{\rm CIP}^2P_{\mathcal{R}\mathcal{R}}$ (where $P_{\mathcal{R}\mathcal{R}}$ is the adiabatic power spectrum) we find, for a number of fiducial galaxy samples in a simplified forecast setup, that constraints on $A_{\rm CIP}$, relative to those on $f_{\rm nl}$, of order $\sigma_{A_{\rm CIP}}/\sigma_{f_{\rm nl}} \approx 1-2$ are achievable for CIPs correlated with adiabatic perturbations, and $\sigma_{A_{\rm CIP}}/\sigma_{f_{\rm nl}} \approx 5$ for the uncorrelated case. These values are independent of survey volume, and suggest that current galaxy data are already able to improve significantly on the tightest existing constraints on CIPs from the CMB. Future galaxy surveys that aim to achieve $\sigma_{f_{\rm nl}} \sim 1$ have the potential to place even stronger bounds on CIPs.

The large-scale clustering of matter is impacted by baryonic physics, particularly AGN feedback. Modelling or mitigating this impact will be essential for making full use of upcoming measurements of cosmic shear and other large-scale structure probes. We study baryonic effects on the matter bispectrum, using measurements from a selection of state-of-the-art hydrodynamical simulations: IllustrisTNG, Illustris, EAGLE, and BAHAMAS. We identify a low-redshift enhancement of the bispectrum, peaking at $k\sim 3h\,{\rm Mpc}^{-1}$, that is present in several simulations, and discuss how it can be associated to the evolving nature of AGN feedback at late times. This enhancement does not appear in the matter power spectrum, and therefore represents a new source of degeneracy breaking between two- and three-point statistics. In addition, we provide physical interpretations for other aspects of these measurements, and make initial comparisons to predictions from perturbation theory, empirical fitting formulas, and the response function formalism. We publicly release our measurements (including estimates of their uncertainty due to sample variance) and bispectrum measurement code as resources for the community.

We introduce The Novel Probes Project, an initiative to advance the field of astrophysical tests of the dark sector by creating a forum that connects observers and theorists. This review focuses on tests of gravity and is intended to be of use primarily to observers, but also to theorists with interest in the development of experimental tests. It is twinned with a separate review on tests of dark matter self-interactions (Adhikari et al., in prep.). Our focus is on astrophysical probes of gravity in the weak-field regime, ranging from stars to quasilinear cosmological scales. These are complementary to both strong-field tests and background and linear probes in cosmology. In particular, the nonlinear screening mechanisms that are an integral part of viable modified gravity models lead to characteristic signals specifically on astrophysical scales. The constraining power of these signals is not limited by cosmic variance, but comes with the challenge of building robust theoretical models of the nonlinear dynamics of stars, galaxies, clusters and large scale structure. In this review we lay the groundwork for a thorough exploration of the astrophysical regime with an eye to using the current and next generation of observations for tests of gravity. We begin by setting the scene for how theories beyond General Relativity are expected to behave, focusing primarily on screened fifth forces. We describe the analytic and numerical techniques for exploring the pertinent astrophysical systems, as well as the signatures of modified gravity. With these in hand we present a range of observational tests, and discuss prospects for future measurements and theoretical developments.

We present a general method to compute the nonlinear matter power spectrum for dark energy and modified gravity scenarios with percent-level accuracy. By adopting the halo model and nonlinear perturbation theory, we predict the reaction of a $\Lambda$CDM matter power spectrum to the physics of an extended cosmological parameter space. By comparing our predictions to $N$-body simulations we demonstrate that with no-free parameters we can recover the nonlinear matter power spectrum for a wide range of different $w_0$-$w_a$ dark energy models to better than 1% accuracy out to $k \approx 1 \, h \, {\rm Mpc}^{-1}$. We obtain a similar performance for both DGP and $f(R)$ gravity, with the nonlinear matter power spectrum predicted to better than 3% accuracy over the same range of scales. When including direct measurements of the halo mass function from the simulations, this accuracy improves to 1%. With a single suite of standard $\Lambda$CDM $N$-body simulations, our methodology provides a direct route to constrain a wide range of non-standard extensions to the concordance cosmology in the high signal-to-noise nonlinear regime.

We present a calculation of the angle-averaged squeezed matter bispectrum covariance ${\rm Cov}\left(B_{m}(k_1, k_1', s_1), B_{m}(k_2, k_2', s_2)\right)$, $s_i \ll k_i,k_i'$ ($i=1,2$), that uses matter power spectrum responses to describe the coupling of large- to short-scale modes in the nonlinear regime. The covariance is given by a certain configuration of the 6-point function, which we show is dominated by response-type mode-coupling terms in the squeezed bispectrum limit. The terms that are not captured by responses are small, effectively rendering our calculation complete and predictive for linear $s_1,s_2$ values and any nonlinear values of $k_1,k_1',k_2,k_2'$. Our numerical results show that the squeezed bispectrum super-sample covariance is only a negligible contribution. We also compute the power spectrum-bispectrum cross-covariance using responses. Our derivation for the squeezed matter bispectrum is the starting point to calculate analytical covariances for more realistic galaxy clustering and weak-lensing applications. It can also be used in cross-checks of numerical ensemble estimates of the general bispectrum covariance, given that it is effectively noise-free and complete in the squeezed limit.

Accurate inference of cosmology from weak lensing shear requires an accurate shear power spectrum covariance matrix. Here, we investigate this accuracy requirement and quantify the relative importance of the Gaussian (G), super-sample covariance (SSC) and connected non-Gaussian (cNG) contributions to the covariance. Specifically, we forecast cosmological parameter constraints for future wide-field surveys and study how different covariance matrix components affect parameter bounds. Our main result is that the cNG term represents only a small and potentially negligible contribution to statistical parameter errors: the errors obtained using the G+SSC subset are within $\lesssim 5\%$ of those obtained with the full G+SSC+cNG matrix for a Euclid-like survey. This result also holds for the shear two-point correlation function, variations in survey specifications and for different analytical prescriptions of the cNG term. The cNG term is that which is often tackled using numerically expensive ensembles of survey realizations. Our results suggest however that the accuracy of analytical or approximate numerical methods to compute the cNG term is likely to be sufficient for cosmic shear inference from the next generation of surveys.

We derive the complete super-sample covariance (SSC) of the matter and weak lensing convergence power spectra using the power spectrum response formalism to accurately describe the coupling of super- to sub-survey modes. The SSC term is completely characterized by the survey window function, the nonlinear matter power spectrum and the full first-order nonlinear power spectrum response function, which describes the response to super-survey density and tidal field perturbations. Generalized separate universe simulations can efficiently measure these responses in the nonlinear regime of structure formation, which is necessary for lensing applications. We derive the lensing SSC formulae for two cases: one under the Limber and flat-sky approximations, and a more general one that goes beyond the Limber approximation in the super-survey mode and is valid for curved sky applications. Quantitatively, we find that for sky fractions $f_{\rm sky} \approx 0.3$ and a single source redshift at $z_S=1$, the use of the flat-sky and Limber approximation underestimates the total SSC contribution by $\approx 10\%$. The contribution from super-survey tidal fields to the lensing SSC, which has not been included in cosmological analyses so far, is shown to represent about $5\%$ of the total lensing covariance on multipoles $\ell_1,\ell_2 \gtrsim 300$. The SSC is the dominant off-diagonal contribution to the total lensing covariance, making it appropriate to include these tidal terms and beyond flat-sky/Limber corrections in cosmic shear analyses.

We compare Einstein-Boltzmann solvers that include modifications to General Relativity and find that, for a wide range of models and parameters, they agree to a high level of precision. We look at three general purpose codes that primarily model general scalar-tensor theories, three codes that model Jordan-Brans-Dicke (JBD) gravity, a code that models f(R) gravity, a code that models covariant Galileons, a code that models Ho\v{r}ava-Lifschitz gravity and two codes that model non-local models of gravity. Comparing predictions of the angular power spectrum of the cosmic microwave background and the power spectrum of dark matter for a suite of different models, we find agreement at the sub-percent level. This means that this suite of Einstein-Boltzmann solvers is now sufficiently accurate for precision constraints on cosmological and gravitational parameters.

We present a calculation of the matter power spectrum covariance matrix ${\rm Cov}(\bf{k}_1,\bf{k}_2)$ that uses power spectrum responses to accurately describe the coupling between large- and small-scale modes beyond the perturbative regime. These response functions can be measured with (small-volume) N-body simulations, which is why the response contributions to the covariance remain valid and predictive at all orders in perturbation theory. A novel and key step presented here is the use of responses to compute loop contributions with soft loop momenta, which extends the application of the response approach beyond that of previously considered squeezed $n$-point functions. The calculation presented here does not involve any fitting parameters. When including response-type terms up to 1-loop order in perturbation theory, we find that our calculation captures the bulk of the total covariance as estimated from simulations up to values of $k_1, k_2 \sim 1\ h/{\rm Mpc}$. Moreover, the prediction is guaranteed to be accurate whenever the softer mode is sufficiently linear, ${\rm min}\{k_1,k_2\} \lesssim 0.08\ h/{\rm Mpc}$. We identify and discuss straightforward improvements in the context of the response approach, which are expected to further increase the accuracy of the calculation presented here.

Cosmological models with Galileon gravity are an alternative to the standard $\Lambda {\rm CDM}$ paradigm with testable predictions at the level of its self-accelerating solutions for the expansion history, as well as large-scale structure formation. Here, we place constraints on the full parameter space of these models using data from the cosmic microwave background (CMB) (including lensing), baryonic acoustic oscillations (BAO) and the Integrated Sachs-Wolfe (ISW) effect. We pay special attention to the ISW effect for which we use the cross-spectra, $C_\ell^{\rm T g}$, of CMB temperature maps and foreground galaxies from the WISE survey. The sign of $C_\ell^{\rm T g}$ is set by the time evolution of the lensing potential in the redshift range of the galaxy sample: it is positive if the potential decays (like in $\Lambda {\rm CDM}$), negative if it deepens. We constrain three subsets of Galileon gravity separately known as the Cubic, Quartic and Quintic Galileons. The cubic Galileon model predicts a negative $C_\ell^{\rm T g}$ and exhibits a $7.8\sigma$ tension with the data, which effectively rules it out. For the quartic and quintic models the ISW data also rule out a significant region of the parameter space but permit regions where the goodness-of-fit is comparable to $\Lambda {\rm CDM}$. The data prefers a non zero sum of the neutrino masses ($\sum m_\nu\approx 0.5$eV) with $ \sim \! 5\sigma$ significance in these models. The best-fitting models have values of $H_0$ consistent with local determinations, thereby avoiding the tension that exists in $\Lambda {\rm CDM}$. We also identify and discuss a $\sim \! 2\sigma$ tension that Galileon gravity exhibits with recent BAO measurements. Our analysis shows overall that Galileon cosmologies cannot be ruled out by current data but future lensing, BAO and ISW data hold strong potential to do so.

We introduce a rigorous definition of general power-spectrum responses as resummed vertices with two hard and $n$ soft momenta in cosmological perturbation theory. These responses measure the impact of long-wavelength perturbations on the local small-scale power spectrum. The kinematic structure of the responses (i.e., their angular dependence) can be decomposed unambiguously through a "bias" expansion of the local power spectrum, with a fixed number of physical response coefficients, which are only a function of the hard wavenumber $k$. Further, the responses up to $n$-th order completely describe the $(n+2)$-point function in the squeezed limit, i.e. with two hard and $n$ soft modes, which one can use to derive the response coefficients. This generalizes previous results, which relate the angle-averaged squeezed limit to isotropic response coefficients. We derive the complete expression of first- and second-order responses at leading order in perturbation theory, and present extrapolations to nonlinear scales based on simulation measurements of the isotropic response coefficients. As an application, we use these results to predict the non-Gaussian part of the angle-averaged matter power spectrum covariance ${\rm Cov}^{\rm NG}_{\ell = 0}(k_1,k_2)$, in the limit where one of the modes, say $k_2$, is much smaller than the other. Without any free parameters, our model results are in very good agreement with simulations for $k_2 \lesssim 0.06\ h/{\rm Mpc}$, and for any $k_1 \gtrsim 2 k_2$. The well-defined kinematic structure of the power spectrum response also permits a quick evaluation of the angular dependence of the covariance matrix. While we focus on the matter density field, the formalism presented here can be generalized to generic tracers such as galaxies.

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 Newton-Gauss-Seidel 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 Millennium-resolution 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, large-sized 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 study the imprints that theories of gravity beyond GR can leave on the lensing signal around line of sight directions that are predominantly halo-underdense (called troughs) and halo-overdense. To carry out our investigations, we consider the normal branch of DGP gravity, as well as a phenomenological variant thereof that directly modifies the lensing potential. The predictions of these models are obtained with N-body simulation and ray-tracing methods using the ECOSMOG and Ray-Ramses codes. We analyse the stacked lensing convergence profiles around the underdense and overdense lines of sight, which exhibit, respectively, a suppression and a boost w.r.t. the mean in the field of view. The modifications to gravity in these models strengthen the signal w.r.t. $\Lambda{\rm CDM}$ in a scale-independent way. We find that the size of this effect is the same for both underdense and overdense lines of sight, which implies that the density field along the overdense directions on the sky is not sufficiently evolved to trigger the suppression effects of the screening mechanism. These results are robust to variations in the minimum halo mass and redshift ranges used to identify the lines of sight, as well as to different line of sight aperture sizes and criteria for their underdensity and overdensity thresholds.

We present high signal-to-noise galaxy-galaxy lensing measurements of the BOSS CMASS sample using 250 square degrees of weak lensing data from CFHTLenS and CS82. We compare this signal with predictions from mock catalogs trained to match observables including the stellar mass function and the projected and two dimensional clustering of CMASS. We show that the clustering of CMASS, together with standard models of the galaxy-halo connection, robustly predicts a lensing signal that is 20-40% larger than observed. Detailed tests show that our results are robust to a variety of systematic effects. Lowering the value of $S_{\rm 8}=\sigma_{\rm 8} \sqrt{\Omega_{\rm m}/0.3}$ compared to Planck2015 reconciles the lensing with clustering. However, given the scale of our measurement ($r<10$ $h^{-1}$ Mpc), other effects may also be at play and need to be taken into consideration. We explore the impact of baryon physics, assembly bias, massive neutrinos, and modifications to general relativity on $\Delta\Sigma$ and show that several of these effects may be non-negligible given the precision of our measurement. Disentangling cosmological effects from the details of the galaxy-halo connection, the effects of baryons, and massive neutrinos, is the next challenge facing joint lensing and clustering analyses. This is especially true in the context of large galaxy samples from Baryon Acoustic Oscillation surveys with precise measurements but complex selection functions.

We perform a validation of estimates of the growth rate of structure, described by the parameter combination $f\sigma_8$, in modified gravity cosmologies. We consider an analysis pipeline based on the redshift-space distortion modelling of the clustering wedges statistic of the galaxy correlation function and apply it to mock catalogues of $\Lambda{\rm CDM}$ and the normal branch of DGP cosmologies. We employ a halo occupation distribution approach to construct our mocks, which we ensure resemble the CMASS sample from BOSS in terms of the total galaxy number density and large scale amplitude of the power spectrum monopole. We show that the clustering wedges model successfully recovers the true growth rate difference between DGP and $\Lambda{\rm CDM}$, even for cases with over 40\% enhancement in $f\sigma_8$ compared to $\Lambda{\rm CDM}$. The unbiased performance of the clustering wedges model allows us to use the growth rate values estimated from the BOSS DR12 data to constrain the cross-over scale $r_c$ of DGP gravity to $\left[r_cH_0\right]^{-1} < 0.97$ ($2\sigma$) or $r_c > 3090\ {\rm Mpc}/h$, cutting into the interesting region of parameter space with $r_c \sim H_0^{-1}$ using constraints from the growth of structure alone.

We present a ray tracing code to compute integrated cosmological observables on the fly in AMR N-body simulations. Unlike conventional ray tracing techniques, our code takes full advantage of the time and spatial resolution attained by the N-body simulation by computing the integrals along the line of sight on a cell-by-cell basis through the AMR simulation grid. Moroever, since it runs on the fly in the N-body run, our code can produce maps of the desired observables without storing large (or any) amounts of data for post-processing. We implemented our routines in the RAMSES N-body code and tested the implementation using an example of weak lensing simulation. We analyse basic statistics of lensing convergence maps and find good agreement with semi-analytical methods. The ray tracing methodology presented here can be used in several cosmological analysis such as Sunyaev-Zel'dovich and integrated Sachs-Wolfe effect studies as well as modified gravity. Our code can also be used in cross-checks of the more conventional methods, which can be important in tests of theory systematics in preparation for upcoming large scale structure surveys.

We introduce and demonstrate the power of a method to speed up current iterative techniques for N-body 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 high-density 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.

Self-consistent ${\it N}$-body simulations of modified gravity models are a key ingredient to obtain rigorous constraints on deviations from General Relativity using large-scale 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 5-10~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 quasi-static limit, which is employed in most modified gravity ${\it N}$-body codes, for the Symmetron model for which the most significant non-static 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 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 N-body 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.

We show that Solar System tests can place very strong constraints on K-mouflage models of gravity, which are coupled scalar field models with nontrivial kinetic terms that screen the fifth force in regions of large gravitational acceleration. In particular, the bounds on the anomalous perihelion of the Moon imposes stringent restrictions on the K-mouflage Lagrangian density, which can be met when the contributions of higher-order operators in the static regime are sufficiently small. The bound on the rate of change of the gravitational strength in the Solar System constrains the coupling strength $\beta$ to be smaller than $0.1$. These two bounds impose tighter constraints than the results from the Cassini satellite and Big Bang Nucleosynthesis. Despite the Solar System restrictions, we show that it is possible to construct viable models with interesting cosmological predictions. In particular, relative to $\Lambda$-CDM, such models predict percent-level deviations for the clustering of matter and the number density of dark matter haloes. This makes these models predictive and testable by forthcoming observational missions.

We show that, in the presence of massive neutrinos, the Galileon gravity model provides a very good fit to the current CMB temperature, CMB lensing and BAO data. This model, which we dub ${\nu} \rm{Galileon}$, when assuming its stable attractor background solution, contains the same set of free parameters as $\Lambda\rm{CDM}$, although it leads to different expansion dynamics and nontrivial gravitational interactions. The data provide compelling evidence ($\gtrsim 6\sigma$) for nonzero neutrino masses, with $\Sigma m_\nu \gtrsim 0.4\ {\rm eV}$ at the $2\sigma$ level. Upcoming precision terrestrial measurements of the absolute neutrino mass scale therefore have the potential to test this model. We show that CMB lensing measurements at multipoles $l \lesssim 40$ will be able to discriminate between the ${\nu} \rm{Galileon}$ and $\Lambda\rm{CDM}$ models. Unlike $\Lambda\rm{CDM}$, the ${\nu} \rm{Galileon}$ model is consistent with local determinations of the Hubble parameter. The presence of massive neutrinos lowers the value of $\sigma_8$ substantially, despite of the enhanced gravitational strength on large scales. Unlike $\Lambda\rm{CDM}$, the ${\nu} \rm{Galileon}$ model predicts a negative ISW effect, which is difficult to reconcile with current observational limits.

We determine the concentration-mass relation of 19 X-ray selected galaxy clusters from the CLASH survey in theories of gravity that directly modify the lensing potential. We model the clusters as NFW haloes and fit their lensing signal, in the Cubic Galileon and Nonlocal gravity models, to the lensing convergence profiles of the clusters. We discuss a number of important issues that need to be taken into account, associated with the use of nonparametric and parametric lensing methods, as well as assumptions about the background cosmology. Our results show that the concentration and mass estimates in the modified gravity models are, within the errorbars, the same as in $\Lambda$CDM. This result demonstrates that, for the Nonlocal model, the modifications to gravity are too weak at the cluster redshifts, and for the Galileon model, the screening mechanism is very efficient inside the cluster radius. However, at distances $\sim \left[2-20\right] {\rm Mpc}/h$ from the cluster center, we find that the surrounding force profiles are enhanced by $\sim20-40\%$ in the Cubic Galileon model. This has an impact on dynamical mass estimates, which means that tests of gravity based on comparisons between lensing and dynamical masses can also be applied to the Cubic Galileon model.

We present a comprehensive derivation of linear perturbation equations for different matter species, including photons, baryons, cold dark matter, scalar fields, massless and massive neutrinos, in the presence of a generic conformal coupling. Starting from the Lagrangians, we show how the conformal transformation affects the dynamics. In particular, we discuss how to incorporate consistently the scalar coupling in the equations of the Boltzmann hierarchy for massive neutrinos and the subsequent fluid approximations. We use the recently proposed K-mouflage model as an example to demonstrate the numerical implementation of our linear perturbation equations. K-mouflage is a new mechanism to suppress the fifth force between matter particles induced by the scalar coupling, but in the linear regime the fifth force is unsuppressed and can change the clustering of different matter species in different ways. We show how the CMB, lensing potential and matter power spectra are affected by the fifth force, and find ranges of K-mouflage parameters whose effects could be seen observationally. We also find that the scalar coupling can have the nontrivial effect of shifting the amplitude of the power spectra of the lensing potential and density fluctuations in opposite directions, although both probe the overall clustering of matter. This paper can serve as a reference for those who work on generic coupled scalar field cosmology, or those who are interested in the cosmological behaviour of the K-mouflage model.

We model the linear and nonlinear growth of large scale structure in the Cubic Galileon gravity model, by running a suite of N-body cosmological simulations using the {\tt ECOSMOG} code. Our simulations include the Vainshtein screening effect, which reconciles the Cubic Galileon model with local tests of gravity. In the linear regime, the amplitude of the matter power spectrum increases by $\sim 25%$ with respect to the standard $\Lambda$CDM model today. The modified expansion rate accounts for $\sim 20%$ of this enhancement, while the fifth force is responsible for only $\sim 5%$. This is because the effective unscreened gravitational strength deviates from standard gravity only at late times, even though it can be twice as large today. In the nonlinear regime ($k \gtrsim 0.1 h\rm{Mpc}^{-1}$), the fifth force leads to only a modest increase ($\lesssim 8%$) in the clustering power on all scales due to the very efficient operation of the Vainshtein mechanism. Such a strong effect is typically not seen in other models with the same screening mechanism. The screening also results in the fifth force increasing the number density of halos by less than 10%, on all mass scales. Our results show that the screening does not ruin the validity of linear theory on large scales which anticipates very strong constraints from galaxy clustering data. We also show that, whilst the model gives an excellent match to CMB data on small angular scales ($l \gtrsim 50$), the predicted integrated Sachs-Wolf effect is in tension with Planck/WMAP results.

We study the nonlinear growth of structure in nonlocal gravity models with the aid of N-body simulation and the spherical collapse and halo models. We focus on a model in which the inverse-squared of the d'Alembertian operator acts on the Ricci scalar in the action. For fixed cosmological parameters, this model differs from $\Lambda{\rm CDM}$ by having a lower late-time expansion rate and an enhanced and time-dependent gravitational strength ($\sim 6\%$ larger today). Compared to $\Lambda{\rm CDM}$ today, in the nonlocal model, massive haloes are slightly more abundant (by $\sim 10\%$ at $M \sim 10^{14} M_{\odot}/h$) and concentrated ($\approx 8\%$ enhancement over a range of mass scales), but their linear bias remains almost unchanged. We find that the Sheth-Tormen formalism describes the mass function and halo bias very well, with little need for recalibration of free parameters. The fitting of the halo concentrations is however essential to ensure the good performance of the halo model on small scales. For $k \gtrsim 1 h/{\rm Mpc}$, the amplitude of the nonlinear matter and velocity divergence power spectra exhibits a modest enhancement of $\sim 12\%$ to $15\%$, compared to $\Lambda{\rm CDM}$ today. This suggests that this model might only be distinguishable from $\Lambda{\rm CDM}$ by future observational missions. We point out that the absence of a screening mechanism may lead to tensions with Solar System tests due to local time variations of the gravitational strength, although this is subject to assumptions about the local time evolution of background averaged quantities.

We use the latest CMB data from Planck, together with BAO measurements, to constrain the full parameter space of Galileon gravity. We constrain separately the three main branches of the theory known as the Cubic, Quartic and Quintic models, and find that all yield a very good fit to these data. Unlike in $\Lambda{\rm CDM}$, the Galileon model constraints are compatible with local determinations of the Hubble parameter and predict nonzero neutrino masses at over $5\sigma$ significance. We also identify that the low-$l$ part of the CMB lensing spectrum may be able to distinguish between $\Lambda{\rm CDM}$ and Galileon models. In the Cubic model, the lensing potential deepens at late times on sub-horizon scales, which is at odds with the current observational suggestion of a positive ISW effect. Compared to $\Lambda$CDM, the Quartic and Quintic models predict less ISW power in the low-$l$ region of the CMB temperature spectrum, and as such are slightly preferred by the Planck data. We illustrate that residual local modifications to gravity in the Quartic and Quintic models may render the Cubic model as the only branch of Galileon gravity that passes Solar System tests.

The velocity field of dark matter and galaxies reflects the continued action of gravity throughout cosmic history. We show that the low-order moments of the pairwise velocity distribution, $v_{12}$, are a powerful diagnostic of the laws of gravity on cosmological scales. In particular, the projected line-of-sight galaxy pairwise velocity dispersion, $\sigma_{12}(r)$, is very sensitive to the presence of modified gravity. Using a set of high-resolution N-body simulations we compute the pairwise velocity distribution and its projected line-of-sight dispersion for a class of modified gravity theories: the chameleon \fR gravity and Galileon gravity (cubic and quartic). The velocities of dark matter halos with a wide range of masses would exhibit deviations from General Relativity at the $(5-10)\sigma$ level. We examine strategies for detecting these deviations in galaxy redshift and peculiar velocity surveys. If detected, this signature would be a "smoking gun" for modified gravity.

We investigate the performance of semi-analytical modelling of large-scale structure in Galileon gravity cosmologies using results from N-body simulations. We focus on the Cubic and Quartic Galileon models that provide a reasonable fit to CMB, SNIa and BAO data. We demonstrate that the Sheth-Tormen mass function and linear halo bias can be calibrated to provide a very good fit to our simulation results. We also find that the halo concentration-mass relation is well fitted by a power law. The nonlinear matter power spectrum computed in the halo model approach is found to be inaccurate in the mildly nonlinear regime, but captures reasonably well the effects of the Vainshtein screening mechanism on small scales. In the Cubic model, the screening mechanism hides essentially all of the effects of the fifth force inside haloes. In the case of the Quartic model, the screening mechanism leaves behind residual modifications to gravity, which make the effective gravitational strength time-varying and smaller than the standard value. Compared to normal gravity, this causes a deficiency of massive haloes and leads to a weaker matter clustering on small scales. For both models, we show that there are realistic halo occupation distributions of Luminous Red Galaxies that can match both the observed large-scale clustering amplitude and the number density of these galaxies.

We study spherical collapse in the Quartic and Quintic Covariant Galileon gravity models within the framework of the excursion set formalism. We derive the nonlinear spherically symmetric equations in the quasi-static and weak-field limits, focusing on model parameters that fit current CMB, SNIa and BAO data. We demonstrate that the equations of the Quintic model do not admit physical solutions of the fifth force in high density regions, which prevents the study of structure formation in this model. For the Quartic model, we show that the effective gravitational strength deviates from unity at late times ($z \lesssim 1$), becoming larger if the density is low, but smaller if the density is high. This shows that the Vainshtein mechanism at high densities is not enough to screen all of the modifications of gravity. This makes halos that collapse at $z \lesssim 1$ feel an overall weaker gravity, which suppresses halo formation. However, the matter density in the Quartic model is higher than in standard $\Lambda$CDM, which boosts structure formation and dominates over the effect of the weaker gravity. In the Quartic model there is a significant overabundance of high-mass halos relative to $\Lambda$CDM. Dark matter halos are also less biased than in $\Lambda$CDM, with the difference increasing appreciably with halo mass. However, our results suggest that the bias may not be small enough to fully reconcile the predicted matter power spectrum with LRG clustering data.

We develop a numerical algorithm to solve the high-order nonlinear derivative-coupling equation associated with the quartic Galileon model, and implement it in a modified version of the RAMSES N-body code to study the effect of the Galileon field on the large-scale matter clustering. The algorithm is tested for several matter field configurations with different symmetries, and works very well. This enables us to perform the first simulations for a quartic Galileon model which provides a good fit to the cosmic microwave background (CMB) anisotropy, supernovae and baryonic acoustic oscillations (BAO) data. Our result shows that the Vainshtein mechanism in this model is very efficient in suppressing the spatial variations of the scalar field. However, the time variation of the effective Newtonian constant caused by the curvature coupling of the Galileon field cannot be suppressed by the Vainshtein mechanism. This leads to a significant weakening of the strength of gravity in high-density regions at late times, and therefore a weaker matter clustering on small scales. We also find that without the Vainshtein mechanism the model would have behaved in a completely different way, which shows the crucial role played by nonlinearities in modified gravity theories and the importance of performing self-consistent N-body simulations for these theories.

We present the first constraints on the full parameter space of the Galileon modified gravity model, considering both the cosmological parameters and the coefficients which specify the additional terms in the Lagrangian due to the Galileon field, which we call the Galileon parameters. We use the latest cosmic microwave background measurements, along with distance measurements from supernovae and baryonic acoustic oscillations, performing a Monte Carlo Markov Chain exploration of the 9-dimensional parameter space. The integrated Sachs-Wolfe signal can be very different in Galileon models compared to standard gravity, making it essential to use the full CMB data rather than the CMB distance priors. We demonstrate that meaningful constraints are only possible in the Galileon parameter space after taking advantage of a scaling degeneracy. We find that the Galileon model can fit the WMAP 9-year results better than the standard \Lambda-Cold Dark Matter model, but gives a slightly worse fit overall once lower redshift distance measurements are included. The best-fitting cosmological parameters (e.g. matter density, scalar spectral index, fluctuation amplitude) can differ by more than 2\sigma\ in the Galileon model compared with \Lambda CDM. We highlight other potential constraints of the Galileon model using galaxy clustering and weak lensing measurements.

We study the cosmology of Galileon modified gravity models in the linear perturbation regime. We derive the fully covariant and gauge invariant perturbed field equations using two different methods, which give consistent results, and solve them using a modified version of the {\tt CAMB} code. We find that, in addition to modifying the background expansion history and therefore shifting the positions of the acoustic peaks in the cosmic microwave background (CMB) power spectrum, the Galileon field can cluster strongly from early times, and causes the Weyl gravitational potential to grow, rather than decay, at late times. This leaves clear signatures in the low-$l$ CMB power spectrum through the modified integrated Sachs-Wolfe effect, strongly enhances the linear growth of matter density perturbations and makes distinctive predictions for other cosmological signals such as weak lensing and the power spectrum of density fluctuations. The quasi-static approximation is shown to work quite well from small to the near-horizon scales. We demonstrate that Galileon models display a rich phenomenology due to the large parameter space and the sensitive dependence of the model predictions on the Galileon parameters. Our results show that some Galileon models are already ruled out by present data and that future higher significance galaxy clustering, ISW and lensing measurements will place strong constraints on Galileon gravity.

We derive the Layzer-Irvine equation in the presence of a homogeneous (or quasi-homogeneous) dark energy component with an arbitrary equation of state. We extend the Layzer-Irvine equation to homogeneous and isotropic universes with an arbitrary number of dimensions and obtain the corresponding virial relation for sufficiently relaxed objects. We find analogous equations describing the dynamics of cosmic string loops and other p-branes of arbitrary dimensionality, discussing the corresponding relativistic and non-relativistic limits. Finally, we generalize the Layzer-Irvine equation to account for a non-minimal interaction between dark matter and dark energy, discussing its practical use as a signature of such an interaction.

In this paper we compare the performances of the chi-square and median likelihood analysis in the determination of cosmological constraints using type Ia supernovae data. We perform a statistical analysis using the 307 supernovae of the Union 2 compilation of the Supernova Cosmology Project and find that the chi-square statistical analysis yields tighter cosmological constraints than the median statistic if only supernovae data is taken into account. We also show that when additional measurements from the Cosmic Microwave Background and Baryonic Acoustic Oscillations are considered, the combined cosmological constraints are not strongly dependent on whether one applies the chi-square statistic or the median statistic to the supernovae data. This indicates that, when complementary information from other cosmological probes is taken into account, the performances of the chi-square and median statistics are very similar, demonstrating the robustness of the statistical analysis.

In this paper we investigate possible solutions to the coincidence problem in flat phantom dark energy models with a constant dark energy equation of state and quintessence models with a linear scalar field potential. These models are representative of a broader class of cosmological scenarios in which the universe has a finite lifetime. We show that, in the absence of anthropic constraints, including a prior probability for the models inversely proportional to the total lifetime of the universe excludes models very close to the $\Lambda {\rm CDM}$ model. This relates a cosmological solution to the coincidence problem with a dynamical dark energy component having an equation of state parameter not too close to -1 at the present time. We further show, that anthropic constraints, if they are sufficiently stringent, may solve the coincidence problem without the need for dynamical dark energy.