There is a well known degeneracy between the enhancement of the growth of large-scale structure produced by modified gravity models and the suppression due to the free-streaming of massive neutrinos at late times. This makes the matter power-spectrum alone a poor probe to distinguish between modified gravity and the concordance $\Lambda$CDM model when neutrino masses are not strongly constrained. In this work, we investigate the potential of using redshift-space distortions (RSD) to break this degeneracy when the modification to gravity is scale-dependent in the form of Hu-Sawicki $f(R)$. We find that if the linear growth rate can be recovered from the RSD signal, the degeneracy can be broken at the level of the dark matter field. However, this requires accurate modelling of the non-linearities in the RSD signal, and we here present an extension of the standard perturbation theory-based model for non-linear RSD that includes both Hu-Sawicki $f(R)$ modified gravity and massive neutrinos.

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 semi-analytical model to calculate SNIe light curves in non-standard 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 (SNIa) 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 can be independently determined, the SNIa can be treated as a tracker to constrain the possible time variation of $G$ in different redshift ranges. The gravitational-wave (GW) standard sirens, caused by the coalescence of binary neutron stars, provide a model-independent way to measure the distance of GW events, which can be used to determine the luminosity distances of SNIa by interpolation, provided the GW and SNIa samples have similar redshift ranges. We demonstrate that combining the GW observations of third-generation detectors with SNIa data provides a powerful and model-independent 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\%$.

The effect of massive neutrinos on the growth of cold dark matter perturbations acts as a scale-dependent Newton's constant and leads to scale-dependent growth factors just as we often find in models of gravity beyond General Relativity. We show how to compute growth factors for $\Lambda$CDM and general modified gravity cosmologies combined with massive neutrinos in Lagrangian perturbation theory for use in COLA and extensions thereof. We implement this together with the grid-based massive neutrino method of Brandbyge and Hannestad in $\texttt{MG-PICOLA}$ and compare COLA simulations to full $\it N$-body simulations of $\Lambda$CDM and $f(R)$ gravity with massive neutrinos. Our implementation is computationally cheap if the underlying cosmology already has scale-dependent growth factors and it is shown to be able to produce results that match $\it N$-body to percent level accuracy for both the total and CDM matter power-spectra up to $k\lesssim 1 h/$Mpc.

We present a general parallelized and easy-to-use code to perform numerical simulations of structure formation using the COLA (COmoving Lagrangian Acceleration) method for cosmological models that exhibit scale-dependent growth at the level of first and second order Lagrangian perturbation theory. For modified gravity theories we also include screening using a fast approximate method that covers all the main examples of screening mechanisms in the literature. We test the code by comparing it to full simulations of two popular modified gravity models, namely $f(R)$ gravity and nDGP, and find good agreement in the modified gravity boost-factors relative to $\Lambda$CDM even when using a fairly small number of COLA time steps.