Institute of Cosmology and Gravitation University of Portsmouth  Kazuya Koyama
Albert Izard
David Wands 0.25 Doctor of Philosophy Cosmology

# Testing gravity on cosmological scales

Theoretical predictions
with the COLA method

Bartolomeo Fiorini

###### Abstract

The standard model of cosmology based on general relativity needs a dark
energy component in the cosmic inventory to successfully describe the
accelerated expansion of the universe at late time. This dark energy is not
well justified by the standard model of particle physics. An alternative way
to explain the accelerated expansion is to consider cosmological models based
on modified gravity theories. The simplest extensions to general relativity
include an additional scalar field mediating a fifth force. The gravity model
is already well constrained with Solar System experiments, which show no sign
of departure from general relativity in this environment. However, gravity
still needs to be precisely constrained on cosmological scales. Among the
modified gravity theories proposed in the literature, we focus on the ones
which incorporate screening mechanisms, hiding the fifth force in high-density
environments.

One of the main objectives of stage IV galaxy surveys is to constrain gravity
on cosmological scales but to fully take advantage of the constraining power
of galaxy surveys it is important to formulate theoretical predictions in the
non-linear regime of structure formation. This is possible by means of N-body
simulations in combination with models for the galaxy-halo connection.
However, full N-body simulations are computationally very expensive and
modified gravity further increases their computational cost. In general
relativity, approximate simulations method can be used to produce synthetic
galaxy catalogues in place of full N-body simulations at the cost of reduced
accuracy in the deep non-linear regime. One of these approximate methods, the
so-called COLA (COmoving Lagrangian Acceleration) method, has been extended to
simulate modified gravity theories.

In this context, we focus on two modified gravity theories, $f(R)$ and nDGP
(the normal branch of the Dvali-Gabadadze-Porrati model), and develop a
pipeline based on the COLA method to produce synthetic galaxy catalogues in
modified gravity. By performing a comparison of COLA summary statistics with
full N-body results, we validate each step of the pipeline and assess the
accuracy of the COLA method. Our results show that COLA is able to accurately
catch the modified gravity effect on the clustering of galaxies in redshift
space, where traces of modified gravity are present in spite of the tuning of
the free parameters of the galaxy-halo connection model. We then use the mock
galaxy catalogues produced with our pipeline to study the effects of modified
gravity and to validate COLA simulations for additional probes of the large-
scale structures. These include an estimator of the power spectrum orthogonal
to the line of sight $Q_{0}$, the bispectrum and voids. Comparing $Q_{0}$ with
the real space power spectrum we show that the modified gravity signal
contained in the two summary statistics is consistent and that COLA accurately
reproduces the N-body results. In the bispectrum of dark matter instead, we
find that COLA simulations, due to the screening approximation that they use,
lack the modified gravity signal coming from the fifth force non-linearity in
$f(R)$ theory. Nonetheless, we show that this is not a problem for the
monopole of bispectrum of galaxies in redshift space as this is dominated by
non-linearities in the bias model, i.e., the model connecting dark matter with
galaxies. We then look at how $f(R)$ and nDGP theories affect the profiles of
voids, finding more modified gravity signatures in nDGP than in $f(R)$
theories. By applying a linear model for the redshift space distortions in
voids we are able to recover unbiased estimates for the linear growth rate in
all gravity theories, even when modified gravity is not taken into account in
the redshift space distortion model. We also show that COLA results for voids
are consistent with the N-body results.

While much faster than full N-body simulations, COLA simulations are still too
computationally expensive to be directly used for cosmological inference which
requires $\sim 10^{5}$ evaluations of theoretical predictions. Here emulation
techniques come in help, requiring as little as $\sim 100$ theoretical
predictions to create smooth interpolating functions that cover a wide range
of the theory parameter space. To pave the way in this direction, we study the
convergence of COLA simulations for predictions of the matter power spectrum
by increasing the force, time and mass resolutions employed in the
simulations. We find that to achieve convergence it is necessary to increase
the three resolutions accordingly. Then we explore the possibility of
extending cosmological emulators with COLA simulations using the response
function for the matter power spectrum, i.e., the response of the matter power
spectrum to changes in cosmological parameters. By comparing COLA predictions
of the response function with state-of-the-art cosmological emulators we show
that COLA is more accurate in predicting the response function than the power
spectrum itself. Finally, we demonstrate the potential of COLA simulations for
the extension of cosmological emulators to modified gravity theories producing
a suite of simulations in nDGP gravity that we employ to train an emulator for
the modified gravity boost factor, i.e., the ratio of the power spectrum in
modified gravity with that in general relativity.

The thesis is structured as follows:

  * •

in chapter 1 we give a general overview of the main concepts at the base of
this work, including the cosmological model, modified gravity theories, the
large-scale structure of the universe and N-body simulations;

  * •

in chapter 2 we present the pipeline for the efficient production of mock
galaxy catalogues with the COLA method and validate it by performing an
extensive comparison of summary statistics with full N-body simulations in
modified gravity;

  * •

in chapter 3 we test the validity of the COLA method for additional probes of
the large-scale structure, a real space power spectrum estimator, bispectrum
and voids, and investigate if they can help constrain modified gravity
theories;

  * •

in chapter 4 we study the feasibility of using COLA simulations to accurately
extend cosmological emulators of the matter power spectrum to modified gravity
theories and give an explicit example in the case of nDGP theory producing an
actual emulator;

  * •

in chapter 5 we summarise the main results discussed in this thesis, draw the
conclusions and discuss future prospects.

###### keywords:

Cosmology Modified gravity COLA Large scale structure Tests of gravity
University of Portsmouth

See Dedication/Dedication.pdf

Whilst registered as a candidate for the above degree, I have not been
registered for any other research award. The results and conclusions embodied
in this thesis are the work of the named candidate and have not been submitted
for any other academic award.

Approximate word count: 27409

## Dissemination

This thesis is based on following works

  1. [1]

Fiorini, Koyama, Izard, Winther, Wright, Li; Fast generation of mock galaxy
catalogues in modified gravity models with COLA, JCAP 09 (2021) 021;

  2. [2]

Fiorini, Koyama, Izard; Studying large-scale structure probes of modified
gravity with COLA, JCAP 12 (2022) 028;

  3. [3]

Brando, Fiorini, Koyama, Winther; Enabling matter power spectrum emulation in
beyond-$\Lambda$CDM cosmologies with COLA, JCAP 09 (2022) 051.

My academic authorship information can be found using either the ORCID record
0000-0002-0092-4321 or the arXiv public author identifier fiorini_b_1.

###### Acknowledgements.

In case you are reading this acknowledgements to look for your name, I wish to
sincerely thank you because of the many reasons that we both know and which I
promise I will not forget. If, on the other hand, you are just curious about
who helped me get to this point, I am sorry to disappoint you: I am very
fortunate to have a lot of important people in my life and all have had a part
in getting me this far. Some of these people clearly contributed more, others
less. Selecting those who contributed the most would inevitably require
leaving out someone who was nonetheless crucial. In any case, Numerical
computations were done on the Sciama High Performance Compute (HPC) cluster
which is supported by the ICG, SEPNet and the University of Portsmouth.

## Notation

In this work we will use the natural units

$c=\hslash\equiv 1$

so that times and distances will be both of the same dimension of the inverse
of energy.

Our metric signature will be $(-+++)$.

Greek indices will take the values ${0,1,2,3}$. The latin indices $i,j,k$ will
be used for the space dimensions so they will run over ${1,2,3}$. These
statements hold unless other definitions are explicitly made. We will adopt
the Einstein notation for repeated indices: when an index variable appears
twice in a single term, it implies summation of that term over all the values
of the index.

Our Fourier convention will be

$\tilde{f}(\mathbf{k})=\int
d^{3}xf(\mathbf{x})e^{-i\mathbf{k}\cdot\mathbf{x}}\,.$

We will use the subscript index “$0$” to denote the present-day values of
variables, unless otherwise stated.

The Hubble expansion rate will be indicated through the letter $H$ and its
conformal version with $\mathcal{H}$. The dot derivatives $\dot{x}$ will stand
for time derivative while the prime derivative $x^{\prime}$ will be for the
conformal time derivative. The partial derivative with respect to the variable
$x^{\mu}$ will be written $\partial_{\mu}:\frac{\partial}{\partial
x^{\mu}}=\partial_{\mu}$.

## Acronyms

2LPT | second-order Lagrangian Perturbation Theory
---|---
AMR | Adaptive mesh refinement
BAO | Barion Acoustic Oscillations
CCF | Cross-Correlation Function
CMB | Cosmic Microwave Background
COLA | COmoving Lagrangian Acceleration
DE | Dark Energy
DGP | Dvali-Gabadadze-Porrati
DM | Dark Matter
EE2 | Euclid Emulator 2
F5 | Hu-Sawicki $f(R)$ model with $\left|f_{R0}\right|=10^{-5}$
FLRW | Friedmann-Lemaitre-Robertson-Walker
FoF | Friends-of-Friends
GR | General Relativity
HOD | Halo Occupation Distribution
IC | Initial Conditions
LHS | Latin-Hypercube Sampling
LSS | Large Scale Structure
MG | Modified Gravity
N1 | nDGP model with $H_{0}r_{c}=1$
nDGP | normal branch of DGP theory
NFW | Navarro-Frenk-White
PCA | Principal Components Analysis
PM | Particle Mesh
RSD | Redshift-Space Distortions
SHAM | Sub-Halos Abundance Matching
SO | Spherical Over-density

###### Contents

  1. 1 Introduction
    1. 1.1 Background cosmology
    2. 1.2 Modified Gravity
      1. 1.2.1 Screening mechanism
      2. 1.2.2 Hu-Sawicki $f(R)$
      3. 1.2.3 Dvali-Gabadadze-Porrati
    3. 1.3 Large scale structure
      1. 1.3.1 Linear theory
      2. 1.3.2 Lagrangian perturbation theory
      3. 1.3.3 Halos
      4. 1.3.4 The bias model
      5. 1.3.5 Redshift space distortions
      6. 1.3.6 Galaxy surveys
    4. 1.4 N-body simulations
      1. 1.4.1 N-body simulations in modified gravity
  2. 2 Fast production of galaxy mock catalogues in modified gravity
    1. 2.1 Simulations
      1. 2.1.1 COLA method
      2. 2.1.2 MG extension to COLA
      3. 2.1.3 Simulation suites
    2. 2.2 Halos
      1. 2.2.1 rockstar halos
      2. 2.2.2 Friends-of-Friends halos
      3. 2.2.3 Internal structure of halos
    3. 2.3 Galaxy catalogues
      1. 2.3.1 HOD model
      2. 2.3.2 HOD parameters fitting
  3. 3 Large-scale structure probes of modified gravity
    1. 3.1 Power spectrum
      1. 3.1.1 Estimator of the real space power spectrum
      2. 3.1.2 Measurement and results
    2. 3.2 Bispectrum
      1. 3.2.1 Definitions and measurements
      2. 3.2.2 Results for the bispectrum
    3. 3.3 Voids
      1. 3.3.1 Redshift-space distortion model in voids
      2. 3.3.2 Velocity model
      3. 3.3.3 Velocity-dispersion estimator
      4. 3.3.4 Void finding and measurements
      5. 3.3.5 Results
  4. 4 Matter power spectrum emulation
    1. 4.1 Convergence tests
      1. 4.1.1 Emulators
      2. 4.1.2 Comparison with emulators and N-body simulations
      3. 4.1.3 Relative convergence
    2. 4.2 Response function
      1. 4.2.1 Variation of cosmological parameters
      2. 4.2.2 Small variations
      3. 4.2.3 Discussion
    3. 4.3 Modified gravity boost factor
      1. 4.3.1 Accuracy and sensitivity
      2. 4.3.2 Simulations and data sets
      3. 4.3.3 Emulator
  5. 5 Summary and conclusions

###### List of Figures

  1. 1.1  Comparison of linear growth factors as a function of the scale factor $a$ in GR (black dashed line), in N1 (orange line) and in F5 gravity (blue lines). The growth factor in F5 gravity is scale-dependent so we show it at three different scales, $k=0.03$, $0.1$ and $0.3\,h\,{\rm Mpc}^{-1}$ (from light to dark blue lines). The bottom panel shows the fractional difference between the growth factors in MG with that in GR.
  2. 1.2  Comparison of matter power spectra at redshift $z=0$ obtained with the N-body codes listed in the legend. Adapted from Angulo:2020vky.
  3. 2.1  Dark Matter power spectrum comparison between COLA (orange dashed lines) and N-body (blue solid lines) obtained by taking the average over 5 realisations. The shaded regions represent the standard deviation over the 5 realisations. The power spectrum in GR (top panels) and the boost-factors in F5 (bottom left) and N1 (bottom right) show agreement within the variance up to $k\sim 1\,h\,{\rm Mpc}^{-1}$.
  4. (a) GR
  5. (b) GR ratio
  6. (c) F5 boost factor
  7. (d) N1 boost factor
  8. 2.2  Dark Matter velocity divergence power spectrum comparison between COLA (orange dashed lines) and N-body (blue solid lines) obtained by taking the average over 5 realisations. The shaded regions represent the standard deviation over the 5 realisations. The power spectrum in GR (top panels) is in agreement within $3\%$ up to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$ where COLA starts to lose power. The boost-factor in F5 (bottom left) and N1 (bottom right) show better than $5\%$ accuracy even deeper in the non-linear regime, up to $k\sim 2\,h\,{\rm Mpc}^{-1}$.
  9. (a) GR
  10. (b) GR ratio
  11. (c) F5 boost factor
  12. (d) N1 boost factor
  13. 2.3 rockstar halo mass function comparison in GR. The hmf in COLA and N-body (left panel) and their ratio (right panel) show $\sim 25\%$ discrepancy for masses greater than $10^{13}\,h^{-1}\,M_{\odot}$.
  14. (a) GR
  15. (b) GR ratio
  16. 2.4 rockstar power spectrum comparison between the power spectra of COLA and N-body in GR. The discrepancy, which is already clear in the left panel, is shown to reach the $10\%$ with a scale-dependent behaviour in the right panel.
  17. (a) GR
  18. (b) GR ratio
  19. 2.5 Mass conversion between FoF and $M_{200c}$ spherical over-density masses. The conversion factor is evaluated for 6 mass values (left panel), corresponding to the SO mass of halos composed of 100, 600, 1000, 3000, 6000 and 10000 DM particles. All the six points lie within $2\%$ from the value $0.75$ which we select as a constant to convert the masses of the FoF halos to SO masses. In the right panel, we use the $N$-body simulations in GR to compare the FoF hmf resulting after this conversion with the hmf of rockstar catalogues. The ratio between the two (right panel) shows a better than $3\%$ agreement for masses above $10^{13}\,h^{-1}\,M_{\odot}$.
  20. (a) Mass conversion factor
  21. (b) rockstar comparison
  22. 2.6 Halo mass function comparison between COLA (orange dashed lines) and N-body (blue solid lines) obtained by taking the average over 5 realisations. The shaded regions represent the standard deviation over the 5 realisations. The hmf in GR (top panels) show a $\sim 5\%$ accuracy for masses above $10^{13}\,h^{-1}\,M_{\odot}$. The boost factor in N1 (bottom right) in COLA agree with that in N-body within the $5\%$ at all masses, while the boost factor in F5 (bottom left) features an agreement within the $10\%$.
  23. (a) GR
  24. (b) GR ratio
  25. (c) F5 boost factor
  26. (d) N1 boost factor
  27. 2.7 Halo power spectrum comparison between COLA (orange dashed lines) and N-body (blue solid lines) obtained by taking the average over 5 realisations. The shaded regions represent the standard deviation over the 5 realisations. The halo power spectrum in GR (top panels) and the boost factors in F5 and N1 (bottom panels) show agreement within the variance up to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$.
  28. (a) GR
  29. (b) GR ratio
  30. (c) F5 boost factor
  31. (d) N1 boost factor
  32. 2.8 Fit of F5 boost factors for concentration parameter (left panel) and velocity dispersion (right panel). The values measured from rockstar catalogues (orange dots) are used to fit the functions in Eq. (2.5) and Eq. (2.6). Only halos with masses above $10^{13}\,h^{-1}\,M_{\odot}$ are used to perform the fits.
  33. (a) Concentration
  34. (b) Velocity dispersion
  35. 2.9 Halo mass (top panels) and velocity distributions (bottom panels) of central (left panels) and satellite (right panels) galaxies in GR. The velocity distribution shows the number of galaxies in each velocity bin, which corresponds to peculiar velocities for centrals and velocities with respect to the host halo rest frame for satellites. The distributions of galaxies in COLA before (orange dot-dashed lines) and after (blue dashed lines) the fit compared with the distributions of N-body in GR (black solid lines) show that the fit is mostly driven by central galaxies.
  36. (a) Mass Distribution Centrals
  37. (b) Mass Distribution Satellites
  38. (c) Velocity Distribution Centrals
  39. (d) Velocity Distribution Satellites
  40. 2.10  Same as Figure 2.9 but for all the gravity models (GR, F5 and N1 in blue, orange and green respectively) after fitting the HOD parameters, which are shown in Table 2.3.
  41. (a) Mass Distribution Centrals
  42. (b) Mass Distribution Satellites
  43. (c) Velocity Distribution Centrals
  44. (d) Velocity Distribution Satellites
  45. 2.11 Comparison of the monopole of the galaxy power spectrum in redshift space between the COLA (in orange) and N-body (in blue) catalogues. The power spectrum in GR (top panels) and the boost factors in F5 and N1 (bottom panels) show an agreement within the variance between COLA and N-body.
  46. (a) GR
  47. (b) GR ratio
  48. (c) F5 boost factor
  49. (d) N1 boost factor
  50. 2.12 Comparison of the quadrupole of galaxy power spectrum in redshift space between COLA (in orange) and N-body (in blue) catalogues. The power spectrum in GR (top panels) show a $5\%$ agreement up to $k\sim 0.17\,h\,{\rm Mpc}^{-1}$. The boost factors in F5 and N1 (bottom panels) show show a better than $5\%$ agreement up to $k\sim 0.2\,h\,{\rm Mpc}^{-1}$, despite $\sim 10\%$ deviations appear above $k\sim 0.17\,h\,{\rm Mpc}^{-1}$ in GR.
  51. (a) GR
  52. (b) GR ratio
  53. (c) F5 boost factor
  54. (d) N1 boost factor
  55. 2.13 Value of $\chi^{2}$ as defined in Eq. (2.9) in COLA (orange dashed lines) and N-body (blue solid lines) evaluated using the signal in GR for the corresponding simulation technique as the target function, as a function of the small scales cut-off $k_{\rm max}$. The $\chi^{2}$ in F5 (left panel) shows a strong scale dependence, with small values on linear scales and larger values on mildly non-linear scales, highlighting the importance of mildly non-linear scales in constraining MG parameters. The $\chi^{2}$ in N1 (right panel) shows little scale dependence, reflecting the features of Vainshtein screening. The results in COLA and N-body are in good agreement.
  56. (a) F5
  57. (b) N1
  58. 3.1  Comparison of the power spectrum orthogonal to the line-of-sight $Q_{0}$ (dots and error-bars) and the real space power spectrum $P^{r}$ (lines and shaded regions) for halos in COLA (in orange) and N-body simulations (in blue). Top left: Full signal in GR. Top right: Ratio of the GR signal between COLA and N-body simulations. Bottom: Boost factors, the ratio of $Q_{0}$ in F5 and N1 to that in GR.
  59. 3.2 Same as figure 3.1 but for galaxies
  60. 3.3  Comparison of bispectrum (left column) and reduced bispectrum (right column) of the matter distribution in COLA (orange crosses) and N-body (blue dots) simulations. The configurations with $k_{1}\geq k_{2}\geq k_{3}$ are displayed in ascending order of the values of the wavenumbers $k_{1}$, $k_{2}$, and $k_{3}$ respectively. The vertical lines denote the value of $k_{1}$ for the configurations immediately to the right of each line. Top: Full signal in GR. Middle and Bottom: Boost factors in F5 and N1 respectively. In each panel, the bottom sub-panel shows the ratio between the COLA and N-body signals displayed in the top sub-panel.
  61. 3.4  Same as figure 3.3 but for the monopole of bispectrum and reduced bispectrum of galaxies in redshift-space.
  62. 3.5 Comparison of the configuration dependence of the F5 (top row) and N1 (bottom row) boost factors of bispectrum (on the left) and reduced-bispectrum (on the right) of DM in N-body (first and third columns) and COLA simulations (second and fourth columns). The colour bars show the amplitude of the boost factors with blue (red) denoting stronger (weaker) signal in MG with respect to GR. In each panel: The top right, top left and bottom corners of the triangle correspond to the equilateral, squeezed and folded configurations respectively. The squeezed configuration is missing from the figure due to the bin’s width.
  63. (a) Full Bispectrum
  64. (b) Reduced Bispectrum
  65. 3.6 Same as figure 3.5 but for the monopole of the bispectrum and reduced bispectrum of galaxies in redshift-space.
  66. (a) Full Bispectrum
  67. (b) Reduced Bispectrum
  68. 3.7  Comparison of bispectrum (left column) and reduced bispectrum (right column) of halos in COLA (orange crosses) and N-body (blue dots) simulations. The configurations with $k_{1}\geq k_{2}\geq k_{3}$ are displayed in ascending order of the values of the wavenumbers $k_{1}$, $k_{2}$, and $k_{3}$ respectively. The vertical lines denote the value of $k_{1}$ for the configurations immediately to the right of each line. Top: Full signal in GR. Middle and Bottom: Boost factors in F5 and N1 respectively. In each panel, the bottom sub-panel shows the ratio between the COLA and N-body signals displayed in top sub-panel.
  69. 3.8  Same as figure 3.7 but for galaxies in real space.
  70. 3.9 Comparison of the configuration dependence of the F5 (top row) and N1 (bottom row) boost factors of bispectrum (on the left) and reduced-bispectrum (on the right) of halos in N-body (first and third columns) and COLA simulations (second and fourth columns). The colour-bars show the amplitude of the boost factors with blue (red) denoting stronger (weaker) signal in MG with respect to GR. In each panel: The top right, top left and bottom corners of the triangle correspond to the equilateral, squeezed and folded configurations respectively. The squeezed configuration is missing from the figure due to the bin’s width.
  71. (a) Full Bispectrum
  72. (b) Reduced Bispectrum
  73. 3.10 Same as figure 3.9 but for galaxies.
  74. (a) Full Bispectrum
  75. (b) Reduced Bispectrum
  76. 3.11 Comparison of the void size function in COLA (in orange) and N-body simulations (in blue). Top left: Full signal in GR. The vertical lines denote the minimum and average void radii after applying the median cut. Top right: Ratio of the GR signal between COLA and N-body simulations. Bottom: Boost factors in F5 and N1. The panels involving ratios are shown in the range $[20,70]\,h^{-1}\,{\rm Mpc}$ as the signal-to-noise ratio deteriorates outside of this range.
  77. 3.12 Comparison of the integrated matter density profile of voids in COLA (in orange) and N-body simulations (in blue). Top left: Full signal in GR. Top right: Ratio of the GR signal between COLA and N-body simulations. Bottom: Boost factors in F5 and N1. The panels involving ratios are shown up to $r\simeq 60\,h^{-1}\,{\rm Mpc}$ as the signal-to-noise ratio deteriorates for larger separations.
  78. 3.13  Comparison of the galaxy radial-velocity profile around voids in COLA (in orange) and N-body simulations (in blue). Top left: Full signal in GR. Top right: Ratio of the GR signal between COLA and N-body simulations. Bottom: Boost factors in F5 and N1. The panels involving ratios are displayed in the range $[20,70]\,h^{-1}\,{\rm Mpc}$ as the signal-to-noise ratio deteriorates outside of the range.
  79. 3.14  Comparison of the galaxy velocity-dispersion profile around voids in COLA (in orange) and N-body simulations (in blue). Top left: Full signal in GR. Top right: Ratio of the GR signal between COLA and N-body simulations. Bottom: Boost factors in F5 and N1. The panels involving ratios are displayed down to $15\,h^{-1}\,{\rm Mpc}$ as the signal-to-noise ratio deteriorates for smaller separations.
  80. 3.15  Comparison of the galaxy-void real-space CCF in COLA (in orange) and N-body simulations (in blue). Top left: Full signal in GR. Top right: Difference of the GR signal between COLA and N-body simulations. Bottom: Difference between the CCF in MG and GR for F5 and N1 gravity. The comparison uses differences instead of usual ratios since the signal crosses 0.
  81. 3.16  Comparison of the monopole of the CCF of voids and galaxies in redshift-space between N-body (left column) and COLA (right column) simulations. Top panels: The difference between the redshift-space monopole of the CCF and the real-space CCF, $\Delta\xi_{0}^{s}\equiv\xi_{0}^{s}-\xi^{r}$, measured in simulations (black dots and shaded region) is compared with the theory predictions (lines) formulated using the radial-velocity profile from simulations ("RSD model", orange) and using the velocity model of eq. (3.22) ("RSD + vel. model", blue). Middle and bottom panels: The difference between $\Delta\xi_{0}^{s}$ in MG and GR (black dots and shaded regions) is compared to the theory predictions (lines) formulated using different combinations of inputs: integrated-density and velocity-dispersion from GR simulations (blue), integrated-density from MG and velocity-dispersion from GR simulations (violet), integrated-density from GR and velocity-dispersion from MG simulations (red), both integrated-density and velocity-dispersion from MG simulations (green). In all cases, we use the real-space CCF from GR simulations in the theory predictions, even in the MG panels to reduce the noise on the MG signals. The velocity model of eq. (3.22) is evaluated with the reference growth rate from linear theory (see table 3.1) in GR and N1 panels. In the case of F5, we use the reference growth rate from the GR linear theory.
  82. 3.17 Same as figure 3.16 but for the quadrupole.
  83. 3.18  Comparison of the galaxy radial-velocity profile in GR N-body simulations with the linear model in eq. (3.22) using the fiducial growth rate $f(z=0.5057)=0.736$, as a function of the distance from the void centre.
  84. 3.19  Correlation of the covariance matrix estimated from jackknife re-sampling in N-body simulations in GR. The black dashed lines are used to distinguish between the monopole (bottom/left) and quadrupole (top/right) terms. The correlation matrices for the other gravity models as well as the ones for N-body simulations appear visually equivalent.
  85. 4.1  Ratio of Bacco and EE2 power spectra for the cosmology defined in eq (4.1).
  86. 4.2  Comparison of PM simulations matter power spectra with Bacco and Arepo power spectra. The high mass-resolution PM run agrees with Bacco up to $k\sim 2h/$Mpc within $1\%$. Arepo does not agree with Bacco but is in better agreement with the low mass-resolution PM run.
  87. 4.3  Ratio of the matter power spectrum of simulations with 512 particles for the force and time resolutions shown in the legend with that of the high-resolution run.
  88. 4.4  Ratio of the matter power spectrum of simulations with 1024 particles for the force and time resolutions shown in the legend with that of the high-resolution run.
  89. 4.5  Ratio of the matter power spectrum of simulations with 2048 particles for the force and time resolutions shown in the legend with that of the high-resolution run.
  90. 4.6 Non-linear and linear response functions for each variation of the cosmological parameters with respect to the reference cosmology. Blue solid lines are computed from COLA simulations, orange dashed lines are obtained from EE2, and dash-dotted green lines are the linear predictions computed using CLASS.
  91. 4.7 Ratio between the non-linear and linear response functions computed from COLA simulations. We show this quantity for five different redshifts: $z=0$ (solid blue), $z=0.5$ (dotted orange), $z=1$ (dashed green), $z=2$ (dash-dotted red) and $z=3$ (solid-squared purple). We use the same convention in all other figures when we show the results at these five redshifts.
  92. 4.8 Ratio between the non-linear response function computed using COLA and the EE2.
  93. 4.9 Ratio between the non-linear response function computed using Bacco and the EE2.
  94. 4.10 Non-linear and linear response functions computed using COLA, EE2 and hi_class, for small variations of the cosmological parameters.
  95. 4.11 The same ratio as in Figure 4.8, but for small variations of the cosmological parameters.
  96. 4.12 The same ratio as in Figure 4.9, but for small variations of the cosmological parameters.
  97. 4.13  Ratio of the nDGP boost factors at redshift $z=1$ between PM simulations and Arepo simulations in 4 nDGP gravity models as indicated in the legend.
  98. 4.14  Ratio of the nDGP boost factors at redshift $z=1$ between PM simulations started at redshift 19 and Arepo simulations in 4 nDGP gravity models as indicated in the legend. The left panel shows the raw results from simulations, the right panel shows the results after applying the linear correction in eq (4.7).
  99. 4.15  Ratio of the nDGP boost factors in low mass-resolution COLA simulations with that in high-resolution PM simulations for two different force resolutions, $N_{\rm mesh}^{1/3}=512$ and $N_{\rm mesh}^{1/3}=2048$. In different panels, we show results for different nDGP gravity models as indicated in the top-left corner of each panel. The black dashed lines show the typical amplitude of the boost factor in the different gravity models for visual reference.
  100. 4.16  Ratio of the nDGP boost factors obtained with $N_{\rm mesh}^{1/3}=1024$ COLA simulations with respect to the same boost factors but obtained with $N_{\rm mesh}^{1/3}=1536$ COLA simulations. The legend describes the redshift value of each line.
  101. 4.17  nDGP boost factor for three values of the parameter $H_{0}r_{c}$. The legend describes the redshift value of each line.
  102. 4.18  Change on the nDGP boost factor for $H_{0}r_{c}=0.2$ due the variation of each cosmological parameter (specified in each panel) with respect to the reference cosmology. The redshift values of each line are as described in the legend.
  103. 4.19  Same as figure 4.18 but for $H_{0}r_{c}=2$
  104. 4.20  Same as figure 4.18 but for $H_{0}r_{c}=20$
  105. 4.21  On the left. Comparison of the two smoothing techniques for the boost factor of nDGP with $H_{0}r_{c}=2$. On the right. Distribution of the normalised parameters in 10 bins between 0 and 1.
  106. 4.22  Comparison between boost factors of the test set (solid lines) and emulated boost factors (dashed lines) for 30 examples randomly chosen in the test set. The residuals between emulated and true boost factors (coloured solid lines) are shown in the bottom panel together with the mean (dashed black line) and variance ($1\sigma$ and $2\sigma$ shaded regions) of the emulation error on the test-set.
  107. 4.23  Comparison between emulator predictions (dashed lines) and Arepo results (solid lines) of the nDGP boost factor at redshift $z=1$ for 4 values of the parameter $H_{0}r_{c}$: 0.5 (in red), 1 (in orange), 2 (in blue) and 5 (in green). The bottom panel shows the residuals.

###### List of Tables

  1. 2.1.3 Simulation suites
  2. 2.3 HOD parameters and number density of galaxies obtained by minimising the objective function in Eq. (2.9) with the simplex search algorithm, except for the N-body GR column where the HOD parameters are the fiducial parameters taken from Manera:2012sc.
  3. 3.1 Results of the optimisation of eq. (3.46) and eq. (3.47) using theory predictions from GR simulations (top) and using theory predictions from MG simulations (bottom). COLA and N-body results are compared side by side in each gravity model (different columns). The additional row at the top of the table shows the linear theory predictions for the linear growth rate at redshift $z=0.5057$ in GR and N1. The linear theory growth rate in F5 is scale-dependent so we omit it from the table.
  4. 4.1 Bacco training set cosmological parameter space.
  5. 4.2 EE2 parameters.
  6. 4.3 Reference parameters.
  7. 4.4 “Large” variation of parameters.
  8. 4.5 Range of cosmological parameters spanned by the nDGP emulator.
  9. 4.6 List of cosmologies used to create the training set for the nDGP emulator.
  10. 4.7 List of cosmologies used to create the test set for the nDGP emulator.

##  Chapter 1 Introduction

> The content of section 1.2 in this chapter is based on the publication
> Fiorini:2021dzs.

###  1.1 Background cosmology

The standard model of cosmology is based on the assumptions that the universe
is homogeneous and isotropic on large-enough scales and that General
Relativity (GR) is the correct description of gravity on all scales of
cosmological interest. With these assumptions, the metric of space-time can be
described using spherical coordinates with the Friedmann-Lemaitre-Robertson-
Walker (FLRW) metric given by the line element

$ds^{2}=-dt^{2}+a(t)^{2}\left(\frac{1}{1-\kappa\chi^{2}}d\chi^{2}+\chi^{2}d\theta^{2}+\chi^{2}\sin^{2}(\theta)d\phi^{2}\right)$
(1.1)

where $a$ is a time-dependent scale parameter and $\kappa$ describes the
spatial curvature of the universe. The coordinate $\chi$ is a comoving
coordinate and it is linked to the physical coordinate $r$ by

$r\equiv a(t)\chi\,.$ (1.2)

This definition is useful to describe the expansion of the universe by means
of the scale factor. The energy-momentum tensor $T_{\mu\nu}$ is related to the
metric of space-time $g_{\mu\nu}$ through the Einstein equations

$R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=8\pi GT_{\mu\nu}\,,$ (1.3)

where $R_{\mu\nu}$ is the Ricci tensor and $R$ is the Ricci scalar, describing
the curvature of space-time. Solving Einstein’s equations with the metric
(1.1) and using the energy-momentum tensor of a perfect fluid

$T_{\mu}^{\nu}=\left(\begin{array}[]{cccc}-\rho&0&0&0\\\ 0&P&0&0\\\ 0&0&P&0\\\
0&0&0&P\end{array}\right)\,,$ (1.4)

where $\rho$ is the energy density and $P$ is the pressure of the fluid, gives
the Friedmann equations

$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}=\frac{8\pi
G}{3}\rho-\frac{\kappa}{a^{2}}\,,$ (1.5)
$\displaystyle\frac{\ddot{a}}{a}=-\frac{4\pi G}{3}(\rho+3P)\,.$ (1.6)

In the standard model of cosmology the content of the universe can be
described as the sum of three perfect fluids with different equations of
state:

  * •

matter, including normal matter (referred to as baryons in cosmology) and cold
Dark Matter (DM), with $P_{\rm m}=0$;

  * •

radiation, including photons and relativistic neutrinos, with $P_{\rm
r}=\frac{1}{3}\rho_{\rm r}$;

  * •

Dark Energy (DE), in the form of a cosmological constant $\Lambda$ entering
the Einstein equation (1.3) on the right hand side, with
$P_{\Lambda}=-\rho_{\Lambda}$ and $\rho_{\Lambda}=\frac{\Lambda}{8\pi G}$.

Due to the assumptions it makes for the dark components (cosmological constant
and cold DM), the standard model of cosmology is also called $\Lambda$CDM
model. Defining the Hubble rate $H\equiv\frac{\dot{a}}{a}$, it is possible to
express the Friedmann equations as

$\left(\frac{H}{H_{0}}\right)^{2}=\Omega_{\rm m}+\Omega_{\rm
r}+\Omega_{\Lambda}+\Omega_{\kappa}$ (1.7)

where $H_{0}$ is the present-day Hubble rate and the dimensionless energy
density $\Omega_{i}\equiv\frac{\rho_{i}}{\rho_{c}}$ with
$\rho_{c}=\frac{3H_{0}^{2}}{8\pi G}$. We have also defined the dimensionless
energy density of curvature $\Omega_{\kappa}$ in the same way, with
$\rho_{\kappa}=-\frac{3\kappa}{8\pi Ga^{2}}$.

The joint constraints of the Cosmic Microwave Background (CMB) measurements
from the Planck mission Planck:2018vyg and the Baryon Acoustic Oscillations
(BAO) measurements from galaxy surveys Beutler:2011hx; Ross:2014qpa;
BOSS:2016wmc are consistent with a null value of $\Omega_{\kappa}$ so in the
following, we will assume a spatially flat universe, $\Omega_{\kappa}=0$.

The supernovae observations show strong evidence for an accelerated expansion
of the universe at present time, consistent with the presence of a DE
component SupernovaSearchTeam:1998fmf; SupernovaCosmologyProject:1998vns as
can be understood from the Friedmann equation (1.6). This DE component, in
principle, could be explained by the vacuum energy of fields in the standard
model of particle physics. However, the amplitude of such vacuum energy would
be $\sim 120$ orders of magnitude larger than what is currently observed
assuming that the ultraviolet cut-off scale is the Planck scale.

With the assumption of a flat universe, the CMB observations constrain the
present-day values of the dimensionless density parameters and Hubble rate to
be Planck:2018vyg

$\begin{array}[]{cc}\Omega_{{\rm m},0}=0.315\pm
0.007\,,&\Omega_{\Lambda,0}=0.685\pm 0.007\,,\\\ \Omega_{{\rm r},0}\simeq
0.0001\,,&H_{0}=67.4\pm 0.5\,{\rm km}\,{\rm Mpc}^{-1}\,{\rm
s}^{-1}\,.\end{array}$ (1.8)

It is possible to take into account small inhomogeneities on top of the
homogeneous background discussed above. These can be described in a GR context
with scalar perturbations of the FLRW metric and the energy-momentum tensor in
the Newtonian gauge. Having previously set $\kappa=0$, the line element reads

$ds^{2}=a^{2}(\tau)\left[-(1+2\Phi)d\tau^{2}+(1-2\Psi)\delta_{ij}dx^{i}dx^{j}\right]$
(1.9)

where $\tau$ is the conformal time defined by $\frac{dt}{d\tau}\equiv a$,
$\mathbf{x}$ is the comoving position in cartesian coordinates and $\Phi$ and
$\Psi$ are scalar perturbations.

Considering small perturbations in the energy density due to the inhomogeneous
distribution of matter

$\rho(\tau,\mathbf{x})=\overline{\rho}(\tau)+\overline{\rho}_{\rm
m}(\tau)\delta_{\rm m}(\tau,\mathbf{x})\,,$ (1.10)

and solving the time-time component of the Einstein equation, we obtain

$\nabla^{2}\Phi=4\pi Ga^{2}\overline{\rho}_{\rm m}\delta_{\rm m}\,,$ (1.11)

which is known as the Poisson equation for the gravitational potential $\Phi$.

###  1.2 Modified Gravity

The $\Lambda$CDM model based on GR has been very successful in reproducing
various cosmological observations. However, GR still lacks a high energy
completion and the $\Lambda$CDM model requires the existence of a highly fine-
tuned cosmological constant to explain the current accelerated expansion of
the Universe. Furthermore, the $\Lambda$CDM model has recently been questioned
for discrepancies between early-universe and late-universe measurements. In
particular, the CMB measurement of $H_{0}$ Planck:2018vyg are in $\sim
5\sigma$ tension with the value inferred from supernovae Riess:2020fzl.
Another tension affecting the $\Lambda$CDM model is the $\sim 3\sigma$ tension
in $S_{8}$ between CMB and stage 3 galaxy surveys Troster:2019ean;
DES:2017qwj; Hildebrandt:2018yau; HSC:2018mrq. These tensions can be
interpreted as hints for new physics beyond the $\Lambda$CDM model that may be
due to the underlying gravity theory Raveri:2019mxg.

These yet-to-be-solved problems have motivated theorists to formulate
alternatives to GR often referred to as modified gravity (MG) theories. A
milestone for the systematic study of gravity theories is represented by
Lovelock’s theorem stating that Einstein’s equations are the only second-order
local equations of motion for a single metric derivable from the covariant
action in four-dimensional spacetime Lovelock:1971yv; Lovelock:1972vz;
Li:2020uaz. This implies that MG theories need to break at least one of the
assumptions of Lovelock’s theorem and this is one of the main approaches
guiding the theoretical efforts in the search for models beyond GR. The MG
theories of interest to cosmology are the ones showing infrared modifications,
which are often associated with an additional force (referred to as the fifth
force), mediated by a scalar field. The knowledge of physically motivated
gravity models is useful for efficient tests of gravity since it reduces the
outcome-space of deviations from GR with respect to model-agnostic approaches.

Solar System tests put very stringent constraints on gravity Will:2014kxa.
Because of this, theories that feature deviations from GR relevant to
cosmology commonly incorporate a screening mechanism Joyce:2014kja;
Koyama:2015vza that let the theory evade these constraints.

In the following, we give an overview of well-known screening mechanisms and
how they are realised in MG theories. Then we focus on two screening
mechanisms of our interest and we introduce two MG theories that naturally
incorporate these screening mechanisms. For a comprehensive review that covers
the topics of this section see Koyama:2018som.

####  1.2.1 Screening mechanism

It is possible to describe the screening mechanisms using the extended Brans-
Dicke action

$\small S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}\left(\psi
R-\frac{\omega(\psi)}{\psi}(\nabla\psi)^{2}+\mathcal{K}\left[(\nabla\psi)^{2},\left(\nabla^{2}\psi\right)\right]-2U(\psi)\right)+S_{m}\left[g_{\mu
v}\right]$ (1.12)

where the scalar field $\psi$ is linearly coupled with the Ricci scalar $R$
and $S_{m}$ is the action of matter. The screening mechanisms can be expressed
in the extended Brans-Dicke formalism for specific functional forms of
$\omega$, $\mathcal{K}$ and $U$:

  * •

The chameleon mechanism Khoury:2003aq; Khoury:2003rn is obtained through the
scalar field potential $U(\psi)$, which gives a density-dependent effective-
mass to the scalar field, thus limiting the range of action of the fifth
force;

  * •

The dilaton Brax:2010gi and symmetron Hinterbichler:2010es mechanisms rely on
the function $\omega(\psi)$, which controls the coupling of the scalar field
to matter. This coupling can be suppressed in high-density regions for
appropriate choices of $\omega(\psi)$;

  * •

The K-mouflage Babichev:2009ee and Vainshtein Vainshtein:1972sx mechanisms use
the non-linear kinetic term
$\mathcal{K}\left[(\nabla\psi)^{2},\left(\nabla^{2}\psi\right)\right]$ to
suppress the coupling with matter. More specifically, they rely on the non-
linearity in the first and second derivatives of the scalar field $\psi$,
respectively.

In general, the Poisson equation (1.11) is modified by the scalar field as

$\nabla^{2}\Phi=4\pi Ga^{2}\rho-\frac{1}{2}\nabla^{2}\psi\,,$ (1.13)

where $\psi$ is the perturbation over the background value of the scalar
field. To close the system of differential equations, the Poisson equation
must be complemented with the scalar field equation (Klein-Gordon) which
depends on the specific screening mechanism under consideration.

The Klein-Gordon equation for the scalar field under the quasi-static
approximation111In the quasi-static approximation, time derivatives of the
scalar field can be neglected compared with spatial derivatives on sub-horizon
scales Noller:2013wca. is, in the chameleon mechanism

$\left(3+2w_{\rm BD}\right)\nabla^{2}\psi+\partial_{\psi}U(\psi)=-8\pi
G\rho\,,$ (1.14)

and in the Vainshtein mechanism

$\left(3+2w_{\rm
BD}\right)\nabla^{2}\psi+\alpha\left[\left(\nabla^{2}\psi\right)^{2}-\left(\nabla_{i}\nabla_{j}\psi\right)^{2}\right]=-8\pi
G\rho\,.$ (1.15)

In the above expressions, $w_{\rm BD}$ and $\alpha$ are model-dependent
quantities. $w_{\rm BD}$ is constrained to be larger than $\approx 40000$ by
solar system experiments in the absence of screening mechanisms Will:2014kxa.
As we will discuss in section 1.4, full N-body simulations in MG (i.e.,
without approximations) solve the Klein-Gordon equations with multi-grid
techniques in order to fully capture the dynamics of the fifth force. However,
screening approximations based on the spherically symmetric solution have been
developed and are currently employed to speed up MG simulations in approximate
simulation methods Winther:2014cia. We discuss these approximations separately
for chameleon and Vainshtein screening in the rest of this section after
performing a frame transformation.

In the action 1.12 the scalar field $\psi$ is coupled with the Ricci scalar
$R$ and matter is minimally coupled with the metric $g_{\mu\nu}$. This is
referred to as the Jordan frame. It is possible to remove the linear coupling
between the scalar field and the Ricci scalar through the conformal
transformation

$g_{\mu\nu}^{E}=\psi g_{\mu\nu}\,.$ (1.16)

The conformally transformed action can be expressed in terms of a field $\phi$
which is minimally coupled with gravity and satisfies

$\left(\frac{d\phi}{d\psi}\right)^{2}=\frac{1}{16\pi
G}\frac{3+2\omega}{\psi^{2}},\quad A(\phi)=\psi^{-1/2},\quad
V(\phi)=\frac{U(\psi)}{8\pi G\psi^{2}}\,.$ (1.17)

This frame takes the name of Einstein frame. In the Einstein frame, the matter
is coupled with gravity through the effective metric
$A^{2}(\phi)g_{\mu\nu}^{E}$ which produces an additional force due to the
scalar field with coupling $\beta/M_{\rm Pl}$, where

$\beta\equiv M_{\mathrm{pl}}\frac{d\ln A}{d\phi}\,.$ (1.18)

In the ideal case of a static, spherically symmetric source of radius $r_{s}$
with density $\rho(r<r_{s})=\rho_{s}$ embedded in a homogeneous background of
density $\rho(r>r_{s})=\rho_{\infty}$, the Klein-Gordon equation for the
Chameleon mechanism in the Einstein frame is

$\frac{d}{dr}\left(r^{2}\frac{d\phi}{dr}\right)=r^{2}\left(V_{,\phi}+\frac{\beta(\phi)\rho_{s}(r)}{M_{\mathrm{Pl}}}\right)$
(1.19)

which has the solution

$\phi(r)=\phi_{\infty}+\frac{\left(\phi_{s}-\phi_{\infty}\right)r_{s}}{r}e^{-m_{\infty}r},\quad
r>r_{s}\,,$ (1.20)

when the gravitational field of the source $\Phi_{N}$ satisfies

$\frac{\Delta
r_{s}}{r_{s}}\equiv\frac{\left|\phi_{\infty}-\phi_{s}\right|}{2\beta_{\infty}M_{\mathrm{Pl}}\Phi_{N}}\ll
1\,$ (1.21)

where $\phi_{\infty}$ is the background value of the scalar field satisfying
$V_{,\phi_{\infty}}+\frac{\beta_{\infty}\rho_{s}(r)}{M_{\mathrm{Pl}}}=0$,
$\beta_{\infty}\equiv\beta(\phi_{\infty})$ and
$m_{\infty}=\left.\frac{d}{d\phi}\left(V_{,\phi}+\frac{\beta(\phi)\rho_{s}(r)}{M_{\mathrm{Pl}}}\right)\right|_{\phi_{\infty}}$
is the effective mass of the scalar field away from the source. The above
expression defines the screening factor $\frac{\Delta r_{s}}{r_{s}}$, which
can be interpreted as the fraction of the radius that the fifth force is able
to penetrate inside the source. Due to this condition, the chameleon mechanism
is also described as thin-shell screening. In the case the gravitational
potential is not strong enough, the source is unscreened and the fifth force
penetrates the source to the centre, $\frac{\Delta r_{s}}{r_{s}}=1$. We can
incorporate this into the theory by re-defining

$\frac{\Delta
r_{s}}{r_{s}}\equiv\operatorname{Min}\left[\frac{\left|\phi_{\infty}-\phi_{c}\right|}{2\beta_{\infty}M_{\mathrm{Pl}}\Phi_{N}},1\right]\,.$
(1.22)

The force perceived by a unit-mass test particle outside the source is then

$\displaystyle F$
$\displaystyle=\frac{GM}{r^{2}}\left[1+2\beta_{\infty}^{2}\left(\frac{\Delta
r_{s}}{r_{s}}\right)\left(1+m_{\infty}r\right)e^{-m_{\infty}r}\right]$ (1.23)
$\displaystyle\simeq\frac{GM}{r^{2}}\left[1+2\beta_{\infty}^{2}\left(\frac{\Delta
r_{s}}{r_{s}}\right)\right]\quad\text{ for }\quad r\ll m_{\infty}^{-1}$
$\displaystyle=F_{\rm N}\left[1+2\beta_{\infty}^{2}\left(\frac{\Delta
r_{s}}{r_{s}}\right)\right]$

where $M$ is the total mass of the source. If we neglect the non-linearity of
the potential $V(\phi)$, the Klein-Gordon equation for the scalar field
becomes

$\nabla^{2}\phi=a^{2}m^{2}(a)\phi+\frac{\beta(a)a^{2}\bar{\rho}_{m}}{M_{\mathrm{Pl}}}\delta_{m}$
(1.24)

where $m(a)\equiv\left.\frac{\mathrm{d}^{2}\left(V+\rho\ln
A\right)}{\mathrm{d}\phi^{2}}\right|_{\phi_{\infty}}$ is the effective mass of
the scalar field around the cosmological background. It is possible to
artificially implement the screening in eq. (1.24) by replacing $\delta_{m}$
with

$\delta_{m}^{\mathrm{eff}}=\delta_{m}\times\operatorname{Min}\left[\frac{\phi(a)}{2\beta(a)M_{\mathrm{Pl}}\left|\Phi_{N}\right|},1\right]$
(1.25)

obtaining the linearised Klein-Gordon equation for the scalar field with
screening approximation

$\nabla^{2}\phi=a^{2}m^{2}(a)\phi+\frac{\beta(a)a^{2}\bar{\rho}_{m}}{M_{\mathrm{Pl}}}\delta_{m}^{\rm
eff}$ (1.26)

In the same spherically symmetric settings as before, the Klein-Gordon
equation for the Vainshtein mechanism in the Einstein frame is given by

$\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\phi}{dr}\right)+\frac{2}{\Lambda_{s}^{3}}\frac{d}{dr}\left(r\left(\frac{d\phi}{dr}\right)^{2}\right)=\frac{\beta\rho_{m}}{M_{\mathrm{Pl}}}$
(1.27)

which can be integrated to obtain

$\frac{d\phi}{rdr}+\frac{2}{\Lambda_{s}^{3}}\left(\frac{d\phi}{rdr}\right)^{2}=2\beta
M_{\mathrm{Pl}}\frac{GM(r)}{r^{3}}$ (1.28)

This can be solved for $F_{\phi}=-\frac{\beta}{M_{\rm Pl}}\nabla_{r}\phi$ to
get

$F_{\phi}=F_{N}\times 2\beta^{2}\times
2\left(\frac{\sqrt{1+\left(r_{V}/r\right)^{3}}-1}{\left(r_{V}/r\right)^{3}}\right)$
(1.29)

where $r_{V}=\frac{1}{\Lambda_{s}}\left(\frac{2\beta M}{\pi
M_{\mathrm{P}1}}\right)^{1/3}$ is the Vainshtein radius $r_{V}$. The ratio
$\left(r_{V}/r\right)^{3}$ can be interpreted as a ratio of densities.
Defining $\rho(r)=\rho_{m}(<r)\equiv\frac{M(r)}{4\pi/3r^{3}}$ and
$\rho_{\mathrm{crit}}=\frac{3\Lambda_{s}^{3}M_{\mathrm{Pl}}}{8\beta}$, the
screening factor in eq. 1.29 can be expressed as

$2\left(\frac{\sqrt{1+\left(r_{V}/r\right)^{3}}-1}{\left(r_{V}/r\right)^{3}}\right)=\frac{2\left(\sqrt{1+\rho(r)/\rho_{\mathrm{crit}}}-1\right)}{\rho(r)/\rho_{\mathrm{crit}}}$
(1.30)

Neglecting the non-linearities, the Klein-Gordon equation for the scalar field
becomes

$\nabla^{2}\phi=\frac{\beta a^{2}\bar{\rho}_{m}}{M_{\mathrm{Pl}}}\delta_{m}$
(1.31)

where we again implement the screening by replacing $\delta_{m}$ with

$\delta_{m}^{\mathrm{eff}}=\delta_{m}\times\frac{2\left(\sqrt{1+\rho(r)/\rho_{\text{crit
}}}-1\right)}{\rho(r)/\rho_{\text{crit }}}$ (1.32)

obtaining the linearised Klein-Gordon equation for the scalar field with
screening approximation

$\nabla^{2}\phi=a^{2}m^{2}(a)\phi+\frac{\beta(a)a^{2}\bar{\rho}_{m}}{M_{\mathrm{Pl}}}\delta_{m}^{\rm
eff}\,.$ (1.33)

Chameleon screening has been tightly constrained by astrophysical tests
studying the morphology of dwarf galaxies in void regions Burrage:2017qrf;
Desmond:2020gzn, at the point that the viable models incorporating the
chameleon mechanisms are almost undistinguishable from GR for astrophysical
and cosmological interest. The Vainshtein mechanism is more difficult to test
with astrophysical tests, due to the highly efficient nature of this screening
mechanism. The strongest constraints come from testing violations of the
strong equivalence principle using supermassive black-holes Bartlett:2020tjd
and from modelling the redshift space distortions in galaxy surveys
Barreira:2016ovx. However, it is still interesting to see whether we can get
comparable constraints from cosmological observations on non-linear scales.

####  1.2.2 Hu-Sawicki $f(R)$

The class of theories where the modifications of gravity can be described by
an additional term to the Einstein-Hilbert action in the form of a function of
the Ricci scalar is known as $f(R)$ theories Starobinsky:1980te:

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

where $\mathcal{L}_{\mathrm{m}}$ is the lagrangian of matter. A sub-set of
$f(R)$ theories, proposed in Hu:2007nk, is known as the Hu-Sawicki $f(R)$
model and is described by

$f(R)=-M^{2}\frac{c_{1}\left(R/M^{2}\right)^{n}}{c_{2}\left(R/M^{2}\right)^{n}+1}\,,$
(1.35)

where $c_{1}$ and $c_{2}$ are dimensionless parameters, $n>0$ is the power law
exponent and $M\equiv H_{0}\sqrt{\Omega_{m,0}}$ is a mass scale defined for
convenience. This model evades the Solar System constraints by means of the
chameleon mechanism Khoury:2003aq; Khoury:2003rn, where the density controls
the shape of the scalar field potential, determining the background value of
the scalar field and its effective mass (large masses in high-density
environments and vice-versa) and the fifth force is consequently screened. In
this work, we focus on the Hu-Sawicki model with $n=1$ and
$\frac{c_{1}}{c_{2}}=6\frac{\Omega_{\Lambda,0}}{\Omega_{m,0}}$ which is
described by the following expression in the high-curvature limit
($M^{2}/R\rightarrow 0$) Hu:2007nk

$f(R)=-2\Lambda+\left|f_{R0}\right|\frac{\bar{R}_{0}^{2}}{R}\,,$ (1.36)

where $\bar{R}_{0}$ is the background curvature today and
$f_{R0}=-n\frac{6\Omega_{\Lambda,0}}{\Omega_{m,0}c_{2}}\left(\frac{\Omega_{m,0}}{3\left(\Omega_{m,0}+4\Omega_{\Lambda,0}\right)}\right)^{n+1}$
is the value of $f_{R}\equiv df/dR$ today which is used as a model parameter
in place of $c_{2}$. For a small $|f_{R0}|$, the background cosmology can be
approximated as the one given by the $\Lambda$CDM model. In the linear regime,
it is described in Fourier space by the modified Poisson equation for the
gravitational potential $\Phi$ Pogosian:2007sw:

$k^{2}\Phi=-4\pi G_{\rm eff}a^{2}\delta\rho,\quad G_{\rm
eff}=G\left(\frac{4+3a^{2}m^{2}/k^{2}}{3+3a^{2}m^{2}/k^{2}}\right),$ (1.37)

where $\delta\rho$ is the fluctuation around the mean energy density and
$m^{2}=\frac{1}{6\left|f_{R0}\right|}\frac{\bar{R}^{3}}{\bar{R}_{0}^{2}}$ is
the effective mass of the scalar field. Depending on the scale and cosmic
time, the effective gravitational constant $G_{\rm eff}$ in Eq. (1.37) has two
limits

$\displaystyle ma\gg k:G_{\rm eff}\rightarrow G\,,$ (1.38) $\displaystyle
ma\ll k:G_{\rm eff}\rightarrow\frac{4}{3}G\,,$ (1.39)

where in the former case the deviations from GR are suppressed while in the
latter gravity is enhanced by the fifth force.

In the non-linear regime, eq. (1.37) for the Poisson equation does not hold
and it is necessary to solve the non-linear scalar field equation with a
potential. The chameleon screening is determined by the gravitational
potential of the object as well as the environment. As we have seen in
subsection 1.2.1, the effect of Chameleon screening can be approximated as
winther15

$G_{\rm
eff}=\left(1+\frac{1}{3}\frac{k^{2}}{k^{2}+a^{2}m^{2}}\text{Min}\left[1,f_{s}\,\left|\frac{\Phi_{\rm
crit}(a)}{\Phi}\right|\right]\right)G,$ (1.40)

where we introduced the empirical parameter $f_{s}$, that we shall call
screening efficiency, to tune the strength of the screening and

$\Phi_{\rm
crit}(a)=\frac{3f_{R0}}{2}\left(\frac{\Omega_{m,0}+4\Omega_{\Lambda,0}}{\Omega_{m,0}a^{-3}+4\Omega_{\Lambda,0}}\right)^{2}\,.$
(1.41)

The quantity $\Phi_{\rm crit}$ depends only on background quantities and
controls (together with $f_{s}$) the threshold value of the gravitational
potential for which the fifth force is screened in high density regions.

####  1.2.3 Dvali-Gabadadze-Porrati

Another way to modify gravity is by defining a theory in a higher-dimensional
space under the assumption that we are living in a 4D brane of the high
dimensional space-time. This class of model is often known as braneworld
gravity and the simplest of these is the Dvali-Gabadadze-Porrati (DGP) model
Dvali:2000hr, which is characterised by the action:

$S=\frac{1}{16\pi G_{(5)}}\int
d^{5}x\sqrt{-\prescript{(5)}{}{g}}\prescript{(5)}{}{R}+\int
d^{4}x\sqrt{-g}\left(\frac{1}{16\pi G}R+\mathcal{L}_{m}\right)\,,$ (1.42)

where the index (5) indicates the quantity is the 5-D generalisation of a 4-D
quantity. The normal branch of this model (nDGP) still needs an additional
dark energy component to explain the late-time acceleration and the tuning of
the cross-over radius $r_{c}\equiv G_{(5)}/2G$ to recover the standard
cosmology at early times.

Although this model is not theoretically well motivated, it is useful as a toy
model to study the effects of its screening mechanism. Indeed focusing on the
normal branch, under quasi-static conditions, the theory incorporates the
Vainsthein mechanism which screens the fifth force in high curvature
environments but allows for deviation from GR in low curvature environments.
We assume that the background cosmology is the same as the $\Lambda$CDM model
by introducing the appropriate dark energy in the model Schmidt:2009sg;
Bag:2018jle. Linearising the equation of motion, the resulting modified
Poisson equation is Koyama:2005kd:

$\nabla^{2}\Phi=4\pi G_{\rm eff}a^{2}\rho\delta,\quad G_{\rm
eff}=G\left(1+\frac{1}{3\beta}\right),$ (1.43)

where $\beta=1-2Hr_{\mathrm{c}}\left(1+\frac{\dot{H}}{3H^{2}}\right)$ and the
overdot denotes the derivative with respect to the physical time. In the case
of the Vainshtein mechanism, the scalar field equation satisfies a non-linear
equation where non-linearity appears in the second derivative of the scalar
field. The Vainshtein screening may be approximated using a density-dependent
effective Newton constant. For the DGP model, the approximation discussed in
subsection 1.2.1 is given by winther15

$G_{\rm
eff}=\left[1+\frac{f_{s}}{3\beta}\left(\frac{2(\sqrt{1+x}-1)}{x}\right)\right]G,\quad
x=\frac{8(r_{c}H_{0})^{2}\Omega_{m,0}}{9\beta^{2}a^{3}}\frac{\rho}{\overline{\rho}},$
(1.44)

where $\rho$ is the average density within a given radius, $\bar{\rho}$ is the
background density and $f_{s}$ is the screening efficiency, an empirical
parameter that we introduce to avoid over-screening at early time by setting
$f_{s}=a^{3}$. Note that this approximation holds only for spherical mass
distribution and can violate the condition that the theory recovers the linear
prediction on large scales Schmidt:2009yj. Also, the screening depends on the
smoothing radius of the density field $\rho$, which needs to be tuned to match
the exact result.

###  1.3 Large scale structure

The cosmic structure that we observe today on very large scales ($\sim$ Mpc)
in galaxy surveys can be understood as the result of the gravitational
evolution of an “initially” almost homogenous energy distribution across the
visible universe. This very homogenous initial state was characterised by
extremely small inhomogeneities ($\sim 10^{-5}$) that we can observe today in
the CMB. These inhomogeneities have been amplified by the gravitational
evolution driven by collisionless cold dark matter. This section aims to
introduce the key elements of the LSS of the universe that will form the
foundations for the analysis of the following chapters. For more on the topics
discussed in this section see Bertschinger:1998tv; Bernardeau:2001qr;
Dodelson:2003ft; Wechsler:2018pic.

####  1.3.1 Linear theory

We consider the distribution of cold matter in a patch of the universe large
enough that its average density can be approximated with the density of the
observable universe for all practical purposes. We focus on a central region
of the patch so that the inhomogeneities of the nearby patches are so far
apart that they have negligible impact on the gravitational field in the
region of interest. We assume that the universe is flat and filled with cold
matter and non-clustering DE so that the density perturbations are entirely
due to perturbations in the matter distribution. Let
$\overline{\rho}(t)=\overline{\rho}_{\rm m}(t)+\overline{\rho}_{\rm\Lambda}$
be the time-dependent density of the observable universe. The density in the
patch can be described by

$\rho(\mathbf{r},t)=\overline{\rho}(t)+\delta\rho_{\rm m}(\mathbf{r},t)\,.$
(1.45)

The gravitational field $\mathbf{\mathcal{F}}$ due to this density
distribution in the Newtonian approximation is given by

$\mathbf{\mathcal{F}}(\mathbf{r},t)=-G\int\rho(\mathbf{s},t)\frac{(\mathbf{r}-\mathbf{s})}{|\mathbf{r}-\mathbf{s}|^{3}}d^{3}s\,.$
(1.46)

Taking the divergence of both sides of the latter equation yields

$\mathbf{\nabla}_{r}\cdot\mathbf{\mathcal{F}}(\mathbf{r},t)=-4\pi
G\rho(\mathbf{r},t)$ (1.47)

where we have made use of the tensor calculus relation
$\mathbf{\nabla}_{r}\cdot\frac{(\mathbf{r}-\mathbf{s})}{|\mathbf{r}-\mathbf{s}|^{3}}=4\pi\delta^{D}(\mathbf{r}-\mathbf{s})$,
with $\delta^{D}$ the Dirac delta distribution. The gravitational field
$\mathbf{\mathcal{F}}$ is irrotational and can therefore be expressed as the
gradient of a scalar field which we denote $\Phi_{\rm N}$, the Newtonian
potential,

$\mathbf{\mathcal{F}}\equiv-\mathbf{\nabla}\Phi_{\rm N}\,.$ (1.48)

Using the definition (1.48) in eq. (1.47) leads to the Poisson equation

$\nabla_{\mathbf{r}}^{2}\Phi_{\rm N}=4\pi G\rho\,,$ (1.49)

which, being a differential equation, requires boundary conditions to be
solved. Since the universe expands according to the Friedmann equations222This
can be equally thought as descending from GR results or the supernovae
observation SupernovaSearchTeam:1998fmf., an effective gravitational field
pulls the matter away from the centre of the reference frame. This
gravitational field can be deduced from the Hubble law $\mathbf{v}_{\rm
bg}=H\mathbf{r}$, by computing the time derivative of the velocity of a
comoving object

$\mathbf{\mathcal{F}}_{\rm
bg}=\frac{d(H\mathbf{r})}{dt}=\frac{d(\mathcal{H}\mathbf{x})}{a\,d\tau}=\frac{\mathcal{H}^{\prime}}{a}\,\mathbf{x}\,.$
(1.50)

To solve the Poisson equation, we impose the boundary condition that the
gravitational field computed from the gravitational potential matches the
background gravitational field $\mathbf{\mathcal{F}}_{\rm bg}$ on average on
the surface of the patch considered. This fixes the shape of the gravitational
potential modulo a constant term that we set to zero. We express the
gravitational potential as the sum of the background value due to the boundary
conditions and a term due to the inhomogeneities

$\Phi_{\rm N}=\Phi_{\rm bg}(t,r)+\phi(t,\mathbf{r})\,$ (1.51)

where $\Phi_{\rm bg}(t)=-\mathcal{H}^{\prime}\,\frac{x^{2}}{2}$. Converting
the Poisson equation in comoving coordinates and using the first Friedmann
equation it is possible to show that

$\nabla_{\mathbf{x}}^{2}\Phi_{\rm bg}=4\pi Ga^{2}\overline{\rho}$ (1.52)

and the Poisson equation becomes

$\nabla_{\mathbf{x}}^{2}\phi=4\pi Ga^{2}\overline{\rho}_{\rm m}\delta_{\rm m}$
(1.53)

which shows that, in the Newtonian description, the gravitational potential
$\phi$ is sourced by the matter density perturbation $\delta_{\rm
m}\equiv\frac{\delta\rho_{\rm m}}{\overline{\rho}_{\rm m}}$.

The conservation of mass implies that the flux of momentum between a volume
$V$ and the surrounding environment through the surface $S$ is the only
responsible for the variation in time of the total matter enclosed in the
volume $V$

$\frac{\partial}{\partial t}\int_{V}\rho=-\oint_{S}\rho\mathbf{u}\cdot dS$
(1.54)

where $\mathbf{u}$ is the velocity field of matter. Applying the Gauss-
Ostrogradsky theorem to the right-hand side, and exploiting the fact that the
resulting equation must hold for every arbitrary volume $V$, gives the
continuity equation

$\frac{\partial\rho}{\partial
t}+\mathbf{\nabla}_{\mathbf{r}}\cdot(\rho\mathbf{u})=0\,.$ (1.55)

The latter can be expressed in comoving coordinates as333The second term in
eq. (1.56) is due to the different slicing of physical and conformal time
coordinates.

$\frac{\partial\rho}{\partial\tau}-\mathcal{H}\mathbf{x}\cdot\mathbf{\nabla}_{\mathbf{x}}\rho+\mathbf{\nabla}_{\mathbf{x}}\cdot(\rho\mathbf{v}+\rho\mathcal{H}\mathbf{x})=0\,.$
(1.56)

where $v$ is the peculiar velocity, and expanding
$\rho=\overline{\rho}(1+\delta)$, gives at first order in the perturbations

$\frac{\partial\delta}{\partial\tau}+\theta=0\,.$ (1.57)

where $\theta=\mathbf{\nabla}_{\mathbf{x}}\cdot\mathbf{v}$ is the velocity-
divergence potential. Applying equivalent reasoning to the conservation of
momentum in the absence of forces leads to

$\left(\partial_{t}+\mathbf{u}\cdot\mathbf{\nabla}_{\mathbf{r}}\right)\mathbf{u}=0\,,$
(1.58)

which becomes the Euler equation,

$\left(\partial_{t}+\mathbf{u}\cdot\mathbf{\nabla}_{\mathbf{r}}\right)\mathbf{u}=-\mathbf{\nabla}_{\mathbf{r}}\Phi\,,$
(1.59)

when gravitational forces are taken into account. The linearised Euler
equation can be expressed in comoving coordinates as

$\partial_{\tau}\mathbf{v}+\mathcal{H}\mathbf{v}=-\mathbf{\nabla}_{\mathbf{x}}\phi$
(1.60)

Continuity, Euler and Poisson equations form a full set of differential
equations uniquely determining the evolution of perturbations. In particular,
taking the comoving divergence of the Euler equation and using the results of
continuity and Poisson equations on its left and right-hand side respectively,
we obtain the second order differential equation in the density contrast

$\delta^{\prime\prime}+\mathcal{H}\delta^{\prime}-\frac{3}{2}\Omega_{\rm
m}H_{0}^{2}\delta=0\,.$ (1.61)

We decompose the density contrast in two factors, one encapsulating its scale
dependence and the other encapsulating its time dependence,
$\delta(\tau,\mathbf{x})=D(\tau)\delta_{\rm ini}(\mathbf{x})$, and substitute
in the eq. (1.61) to obtain the differential equation

$D^{\prime\prime}+\mathcal{H}D^{\prime}-\frac{3}{2}a^{2}\Omega_{\mathrm{m}}H_{0}^{2}D=0\,.$
(1.62)

The quantity $D(\tau)$ is often referred to as the linear growth factor as it
describes the growth of density perturbation on large scales, where linear
theory is accurate. The differential equation (1.62) has both growing and
decaying modes. From the growing one, $D^{(+)}$, we define the linear growth
rate

$f\equiv\frac{\mathrm{d}\ln D^{(+)}}{\mathrm{d}\ln a}\,,$ (1.63)

that will be useful in subsection 1.3.5 and in section 3.3.

This treatment is valid in a flat universe with cold matter and DE, therefore
it is believed to be a good description of our universe during the era of
matter and DE domination. The presence of relativistic species requires a GR
treatment and affects both the gravitational potential and the background
expansion but their effect on the growth function can be neglected.

When MG is taken into account the Poisson equation is modified by the
introduction of a (possibly) scale-dependent Newton constant in Fourier space
as discussed in section 1.2. The growth equation in MG theories characterised
by an effective Newton constant $G_{\rm eff}$ reads

$D^{\prime\prime}(\tau,k)+\mathcal{H}D^{\prime}(\tau,k)-\frac{3}{2}\frac{G_{\rm
eff}(\tau,k)}{G}a^{2}\Omega_{\mathrm{m}}H_{0}^{2}D(\tau,k)=0\,.$ (1.64)

This gives rise to a scale-dependent linear growth rate if $G_{\rm eff}$ is
scale-dependent. Otherwise, the linear growth rate is scale independent as in
the GR case. The solution of the growth factor as a function of the scale
factor can be computed with numerical methods for a fixed $\Omega_{\rm m}$. In
figure 1.1

Fig. 1.1:  Comparison of linear growth factors as a function of the scale
factor $a$ in GR (black dashed line), in N1 (orange line) and in F5 gravity
(blue lines). The growth factor in F5 gravity is scale-dependent so we show it
at three different scales, $k=0.03$, $0.1$ and $0.3\,h\,{\rm Mpc}^{-1}$ (from
light to dark blue lines). The bottom panel shows the fractional difference
between the growth factors in MG with that in GR.

we show a comparison of linear growth factors in GR (black dashed line), in
the nDGP theory with parameter $H_{0}r_{c}=1$ (N1, orange line) and in the Hu-
Sawicky $f(R)$ theory with $n=1$ and $|f_{R0}|=10^{-5}$ (F5, blue lines),
obtained with $\Omega_{\rm m}=0.281$. Since $G_{\rm eff}$ in $f(R)$ theories
is scale dependent, we show the linear growth factor in F5 at 3 different
scales $k=0.03$, $0.1$ and $0.3\,h\,{\rm Mpc}^{-1}$ (from light to dark blue
lines). At present time ($a=1$) the linear growth factor in N1 is $\approx
7\%$ larger than in GR. Also in F5 the linear growth factor is enhanced at
$a=1$ compared to GR, but the enhancement varies with the scale considered,
from $\approx 0.5\%$ at $k=0.03\,h\,{\rm Mpc}^{-1}$ to $\approx 11\%$ at
$k=0.3\,h\,{\rm Mpc}^{-1}$. This can be understood as the result of the
effective mass of the scalar field in F5 that limits the range of interaction
of the fifth force leaving the growth on larger scales almost unaffected.

To derive eq 1.64, we implicitly used the Fourier transform of the density
contrast

$\tilde{\delta}(\mathbf{k},\tau)=\int\delta(\mathbf{x},\tau)e^{-i\mathbf{k}\cdot\mathbf{x}}\,d^{3}x\,,$
(1.65)

and separated again the scale and time dependencies of the density contrast
$\tilde{\delta}(\mathbf{k},\tau)=D(k,\tau)\tilde{\delta}_{\rm
ini}(\mathbf{k})$ The Fourier transform of the density contrast is useful to
define the power spectrum

$\left\langle\tilde{\delta}(\mathbf{k},\tau)\tilde{\delta}\left(\mathbf{k}^{\prime},\tau\right)\right\rangle\equiv(2\pi)^{3}P(k,\tau)\delta_{D}\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\,,$
(1.66)

a summary statistics which contains important information of the density
field. In the regime of applicability of linear theory, it is possible to show
that the power spectrum scales as the square of the growth factor

$P(k,\tau)=D^{2}(k,\tau)P_{\rm ini}(k)\,,$ (1.67)

where we stress again that $D(k,\tau)$ depends on the scale only for theories
where $G_{\rm eff}$ is scale-dependent.

####  1.3.2 Lagrangian perturbation theory

Instead of studying the evolution of the density contrast of given volume
elements using Eulerian coordinates $\mathbf{x}$ as done in the previous sub-
section, it can be convenient to track the motion of fluid elements identified
with Lagrangian coordinates $\mathbf{q}$ using a displacement field
$\mathbf{\xi}$ Bernardeau:2001qr. The two coordinates’ systems are related by
the mapping

$\mathbf{x}=\mathbf{q}+\mathbf{\xi}(\mathbf{q},\tau)\,.$ (1.68)

The equation of motion for the fluid elements in the Newtonian approximation
is given by

$\frac{\mathrm{d}^{2}\mathbf{x}}{\mathrm{d}\tau^{2}}+\mathcal{H}(\tau)\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\tau}=-\mathbf{\nabla}_{x}\phi$
(1.69)

which can be expressed in terms of the displacement field as

$\frac{\mathrm{d}^{2}\mathbf{\xi}}{\mathrm{d}\tau^{2}}+\mathcal{H}(\tau)\frac{\mathrm{d}\mathbf{\xi}}{\mathrm{d}\tau}=-\mathbf{\nabla}_{x}\phi$
(1.70)

or, using the super-comoving time coordinate $\eta$ such that
$d\eta=ad\tau=a^{2}dt$,

$\frac{\mathrm{d}^{2}\mathbf{\xi}}{\mathrm{d}\eta^{2}}=-\mathbf{\nabla}_{x}\phi_{\rm
sc}\,,$ (1.71)

where we introduced the super-comoving gravitational potential $\phi_{\rm
sc}\equiv a^{2}\phi$. Assuming that the mapping (1.68) is bijective, it is
possible to express the displacement field in terms of the density contrast by
means of the mass conservation relation

$\int_{V}d^{3}x\rho_{\mathrm{m}}(\mathbf{x},\eta)=\int_{V}d^{3}q\bar{\rho}_{\mathrm{m}}(\eta)\quad\Rightarrow\quad(1+\delta)\mathrm{J}_{\mathbf{x},\mathbf{q}}=1$
(1.72)

where $\mathrm{J}_{\mathbf{x},\mathbf{q}}\equiv\det(\frac{\partial
x_{i}}{\partial q_{j}})$. Taking the divergence of eq. (1.71) and using the
Poisson equation (1.53) we get

$\frac{\mathrm{d}^{2}\xi_{i,i}}{\mathrm{d}\eta^{2}}-\xi_{j,i}\frac{d^{2}}{d\eta^{2}}\xi_{i,j}=-4\pi
G\bar{\rho}_{\mathrm{m}}\delta\,,$ (1.73)

where $\xi_{i,j}\equiv\partial\xi_{i}/\partial q_{j}$. Expanding both the
density contrast $\delta$ and the displacement field $\mathbf{\xi}$ in a
perturbative series in terms of a small parameter $\epsilon$,

$\displaystyle\mathbf{\xi}=\epsilon\mathbf{\xi}^{(1)}+\epsilon^{2}\mathbf{\xi}^{(2)}+...\,,$
(1.74)
$\displaystyle\delta=\epsilon\delta^{(1)}+\epsilon^{2}\delta^{(2)}+...\,,$
(1.75)

allows us to find the solutions of the set of equations (1.72) and (1.73)
order by order in perturbation theory. This framework is known as Lagrangian
Perturbation Theory (LPT) Bernardeau:2001qr. We assume that the curl of the
displacement field vanishes so that the displacement field can be expressed as
the gradient of a potential field, $\mathbf{\xi}=\mathbf{\nabla}\phi_{\rm
LPT}$. Factoring out the time dependence of the displacement field as
$\mathbf{\xi}^{(i)}(\eta,\mathbf{q})=D_{(i)}(\eta)\mathbf{\xi}^{(i)}(\eta=\eta_{\rm
ini},\mathbf{q})$, the first order solution is given by

$\displaystyle\nabla_{q}^{2}\phi_{\rm LPT}^{(1)}=-D_{1}(\eta)\delta_{\rm
ini}(\mathbf{q})\,,$ (1.76) $\displaystyle\frac{d^{2}D_{1}}{d\eta^{2}}-4\pi
G\bar{\rho}_{\mathrm{m}}{D_{1}}=0\,,$ (1.77)

Truncating the perturbation theory expansion at the first order gives rise to
the so-called Zel’dovich approximation Zeldovich:1969sb. The second-order
Lagrangian perturbation theory (2LPT) solution is given by

$\displaystyle\nabla_{q}^{2}\phi_{\rm
LPT}^{(2)}=\frac{1}{2}D_{2}(\eta)\sum_{i\neq
j}\left(\xi_{i,i}^{(1)}\xi_{j,j}^{(1)}-\xi_{i,j}^{(1)}\xi_{j,i}^{(1)}\right)$
(1.78) $\displaystyle\frac{d^{2}D_{2}}{d\eta^{2}}-4\pi
G\bar{\rho}_{\mathrm{m}}D_{2}=-4\pi G\bar{\rho}_{\mathrm{m}}D_{1}^{2}\,.$
(1.79)

When MG is taken into account, the Poisson equation is in general modified by
a time and scale dependent term as in eq. (1.13). In this case, the LPT
solutions can be found mode by mode in Fourier space in terms of modified
scale-dependent growth factors Valogiannis:2016ane; Winther:2017jof;
Aviles:2017aor.

####  1.3.3 Halos

The growth factor equation (1.62) describes how the initially small
fluctuations of the density field grow over time in the linear regime.
However, as more matter falls towards the overdensities, the motion of DM
becomes non-linear and bound objects of DM start forming. These are called
halos.

To understand halos it is useful to study the ideal case of a spherical top-
hat overdensity in a homogeneous expanding universe. The magnitude of
gravitational energy of the outermost shell due to the overdensity must be
larger than its kinetic energy for the shell to stop expanding, revert its
motion and collapse reaching the virial equilibrium. In a universe dominated
by cold DM, the problem has an analytical (parametric) solution which predicts
that the spherical overdensity collapses to a point when its linear evolution
would hit the threshold of $\delta_{c}\simeq 1.686$. However, instead of
collapsing to a point the matter distribution reaches the virial equilibrium
with an overdensity $\Delta_{\rm vir}\approx 180$. Including the DE component
into the problem affects the results for $\delta_{c}$ and $\Delta_{\rm vir}$
only marginally as the collapsing structures are decoupled from the Hubble
flow when DE becomes relevant. Despite representing an approximation (or
better the monopole) of the aspherical halo collapse, this model catches some
key properties of halos that are confirmed by computer simulations:

  * •

halo collapse is a local process: the matter belonging to the spherically
collapsed halo was originally inside the outer-most shell of the overdensity,
with initial radius $R_{\rm i}=\sqrt[3]{\frac{3M}{4\pi\rho_{\mathrm{m},0}}}$
for a halo of mass $M$,

  * •

the density is (approximately) the same for all virialised objects: in the
spherical collapse model it is determined by the virial theorem,
$\rho=\Delta_{\rm vir}\times\bar{\rho}_{\mathrm{m}}$.

The spherical collapse framework is used as a starting point in the extended
Press-Schechter formalism Press:1973iz; Bond:1990iw, which applies the
excursion-set theory to the Brownian motion of the density contrast as a
function of its variance $\sigma^{2}(R)$ computed in spheres of comoving
radius $R$. For a single Brownian walk, i.e. studying the density centred in a
single point in space, the first crossing of the barrier $\delta_{c}$
determines the current mass of the halo which hosts the element of mass
initially located in the point considered. From a different perspective,
studying the fraction of walks that are above the threshold $\delta_{c}(z)$ in
terms of the variance $\sigma^{2}$, gives the fraction of halos collapsed as a
function of their mass and redshift in the Press-Schechter formalism. This can
be used to estimate the number density of halos $n(M,z)$ at a given mass and
redshift which satisfies

$\frac{dn(M,z)}{dM}=\frac{\rho_{\mathrm{m},0}}{M^{2}}f_{\mathrm{PS}}\left(\frac{\delta_{\mathrm{c}}}{\sigma(M,z)}\right)\left|\frac{d\ln\sigma(M,z)}{d\ln
M}\right|,\quad f_{\mathrm{PS}}(v)=\sqrt{\frac{2}{\pi}}\nu e^{-v^{2}/2}\,,$
(1.80)

where $M$ is the mass enclosed in the sphere of comoving radius $R$,
$M(R)=\frac{4\pi}{3}R^{3}\rho_{{\rm m},0}$. For large halo masses, which
correspond to spheres of large comoving radii,
$\sigma(M,z)\ll\delta_{\mathrm{c}}$ and the abundance of halos is
exponentially suppressed.

This description is useful to gain intuition about how halos form and why
halos of large mass are very rare. However, this is just a toy model and to
have reliable halo statistics it is important to fully take into account the
non-linear dynamics of halo collapse. This can be done by running high
resolution N-body simulations which simulate the dynamics of a large number of
DM particles under the action of gravity forces (see section 1.4 for more
details). Thanks to the analysis of N-body simulations it is possible to
derive fitting formulae for the average radial profile of halos and for the
halo mass function, $n(>M,z)$, i.e., the abundance of halos with a mass larger
than $M$ as a function of $M$.

A practical way of identifying halos in simulations is to identify the largest
spheres of comoving radius $R_{\Delta}$ that contain the same average density

$\frac{M\left(<R_{\Delta}\right)}{4\pi R_{\Delta}^{3}/3}=\Delta_{\rm
m}\times\rho_{\mathrm{m},0}=\Delta_{\rm c}\times\rho_{\mathrm{c}}$ (1.81)

where $\Delta$ is the overdensity parameter that can be defined in terms of
the mean matter density $\rho_{{\rm m},0}$ or the critical matter density
$\rho_{\rm c}$444The critical overdensity parameter is related to the matter
overdensity parameter by $\Delta_{\rm c}=\Omega_{{\rm m},0}\Delta_{\rm m}$..
The halos found in this way are referred to as spherical overdensity halos and
their properties depend on the overdensity parameter $\Delta$.

The profile of spherical overdensity halos has been shown to be approximately
universal and well described by the so-called Navarro-Frenk-White (NFW)
profile Navarro:1995iw

$\rho_{\mathrm{NFW}}(r)=\frac{\rho_{s}}{\left(r/r_{s}\right)\left(1+r/r_{s}\right)^{2}}$
(1.82)

in terms of two parameters, the scale radius $r_{s}$ and the density
$\rho_{s}$. Alternatively, the NFW profile can be expressed in terms of the
concentration parameter $c_{\Delta}\equiv R_{\Delta}/r_{s}$ and the mass

$M_{\Delta}\equiv\int_{0}^{R_{\Delta}}4\pi r^{2}\rho_{\rm
NFW}(r)dr=4\pi\rho_{s}r_{s}^{3}\left[\ln(1+c)-\frac{c}{1+c}\right]\,.$ (1.83)

The velocity of DM particles inside halos can be described as the sum of two
components: a circular velocity component $v_{\rm circ}$ and a radial velocity
component $v_{r}$. Under the assumption that the matter inside the sphere of
radius $R_{\Delta}$ is virialised, the NFW profile can be used to predict the
average velocity profiles of halos. The average modulus of the circular
velocity depends on the mass $M(<r)$ enclosed in the sphere of radius $r$ as

$v_{\rm circ}(r)=\sqrt{\frac{GM(<r)}{r}}$ (1.84)

while the radial-velocity dispersion profile can be expressed as NFW_VelDisp

$\langle v_{r}^{2}(r)\rangle=\frac{1}{\rho_{\text{NFW
}}(r)}\int_{r}^{\infty}\mathrm{d}r\rho_{\text{NFW
}}(r)\frac{\mathrm{d}\Phi(r)}{\mathrm{d}r}\,,\quad\Phi(r)=\frac{-GM(<r)}{r}.$
(1.85)

The halo mass functions found in N-body simulations have been used in
literature to tune analytical fitting functions (see for instance
Jenkins:2000bv; Sheth:2001dp; Warren:2005ey) and more recently to train
emulators (e.g., Heitmann:2015xma; McClintock:2018uyf; Bocquet:2020tes). In
particular, by estimating the abundance of spherical overdensity halos from a
suite of N-body simulations with different resolutions for a wide range of
halo masses, the authors of Tinker:2008ff have produced a fitting formula for
the halo mass function with $5\%$ accuracy in the mass range $10^{11}\leq
M\leq 10^{15}h^{-1}M_{\odot}$, becoming the standard in cluster abundance
analysis of the last decade.

####  1.3.4 The bias model

Let us consider a long wavelength density perturbation
$\delta_{\ell}(\mathbf{x},t)$ in the excursion set theory. The perturbation
$\delta_{\ell}$ adds up to the smaller scales density contrast smoothed on
some radius $R$, $\delta_{R}(\mathbf{x},t)$, increasing the chances of the
density contrast to hit the spherical collapse barrier $\delta_{\rm c}$ in the
regions of space where $\delta_{\ell}(\mathbf{x},t)>0$. From a different point
of view, the perturbation $\delta_{\ell}(\mathbf{x},t)$ can be seen as
affecting the height of the barrier $\delta_{\rm c}$

$\delta_{\mathrm{c}}\longrightarrow\delta_{\mathrm{c}}-\delta_{\ell}^{(1)}(\boldsymbol{x},t)\,.$
(1.86)

where the apex $(1)$ indicates that the density contrasts is evolved to the
time $t$ using the linear growth factor $D(t)$ discussed in section 1.3.1. Now
the expected number density of halos of eq. (1.80), in the presence of this
long wavelength overdensity, becomes

$\left.\frac{dn(M,z)}{dM}\right|_{\delta_{\ell}}=\frac{\rho_{\mathrm{m}}\left(t_{0}\right)}{M^{2}}f_{\mathrm{PS}}\left(\frac{\delta_{\mathrm{c}}-\delta_{\ell}^{(1)}(\mathbf{x},t)}{\sigma(M,z)}\right)\left|\frac{d\ln\sigma(M,z)}{d\ln
M}\right|\,,$ (1.87)

which defines the halo bias $b(M,z)$ due to the long wavelength perturbation
$\delta_{\ell}^{(1)}(\mathbf{x},t)$ in terms of the halo number density
contrast as

$\delta_{\mathrm{h},\ell}(\mathbf{x},t)=\frac{dn/\left.d\ln
M\right|_{\delta_{\ell}}}{dn/\left.d\ln M\right|_{0}}-1\equiv
b(M,z)\delta_{\ell}^{(1)}(\mathbf{x},t)\,.$ (1.88)

At the first order in $\delta_{\ell}^{(1)}(\boldsymbol{x},t)$, in the extended
Press-Schechter formalism, the linear bias is given by

$b_{1}^{\mathrm{PS}}(M,z)=\frac{\delta_{\rm
c}^{2}-\sigma^{2}(M,z)}{\delta_{\mathrm{c}}\sigma^{2}(M,z)}\approx\frac{\delta_{\rm
c}}{\sigma^{2}(M,z)}\,,\quad\text{for}\quad\sigma(M,z)\ll\delta_{\rm c}$
(1.89)

which shows that the bias is larger for larger halo masses (smaller
$\sigma(M,z)$). From a statistical perspective, the halo bias describes how
much more clustered the halos are with respect to the underlying matter
distribution. This concept of bias can be extended to other tracers, such as
galaxies.

In the current understanding of galaxy formation, galaxies form inside halos
Wechsler:2018pic. Hydrodynamical simulations let us study the formation of
galaxies by incorporating hydrodynamical models of subgrid physics in high-
resolution N-body simulations555See section 1.4 for more details on N-body
simulations. However, hydrodynamical simulations are computationally very
expensive, so empirical prescriptions are often used to connect halos found in
DM-only simulations with galaxies. Among these, we focus on abundance matching
and Halo Occupation Distribution (HOD) models.

The abundance matching prescription assumes that the most massive or luminous
galaxies live in the most massive halos, so if $N_{\rm g}$ galaxies are to be
placed in a catalogue of $N_{\rm h}$ halos, halos are ordered by mass and the
$N_{\rm g}$ most massive halos are selected to host a galaxy. This method has
the advantage of being non-parametric, but to provide accurate results the
halos substructure needs to be well resolved Klypin:2013rsa. This is possible
in high-resolution N-body simulations where both halos and sub-halos are
identified and galaxies are assigned to them as central and satellite galaxies
respectively. When abundance matching is applied to both halos and sub-halos
it is often referred to as Sub-Halos Abundance Matching (SHAM)
Kravtsov:2003sg; Vale:2004yt. One difficulty of accurately employing SHAM is
that subhalos and galaxies are stripped of their mass by the host halos in
different ways Nagai:2004ac, hence galaxies should be matched to sub-halos at
the time they are accreted. This requires tracking the evolution of halos with
merger-trees Conroy:2005aq further increasing the computational cost of the
already expensive high-resolution N-body simulations. Conversely, in
simulations with lower resolution, only host halos are resolved
Wechsler:1997fz and the clustering statistics predicted by this technique are
not accurate enough to be employed in modern galaxy surveys Klypin:1997fb.

The HOD is a parametric method to populate halos with galaxies. It gives the
probability, $P(N_{\rm g}\mid M_{\rm h},\dots)$, that a number of galaxies,
$N_{\rm g}$, are found in a halo depending on its properties, normally on its
mass $M_{\rm h}$. The galaxies are often split between central and satellite
galaxies and their probability distributions follow the Bernoulli and Poisson
distributions respectively. The functional form for the satellite galaxies is
studied in high-resolution simulations assuming that the number of satellite
galaxies depends on the halo mass similarly to the number of subhalos, hence a
power law is often assumed Kravtsov:2003sg. The satellites’ positions and
velocity are assigned based on analytical halo profiles whose functional form
is derived from high-resolution simulations, such as NFW or Einasto profile
Navarro:1995iw; Graham:2005xx. On the one hand, HOD models can be applied to
low-resolution simulations since they rely on basic halo properties (e.g., the
halo mass). On the other hand, they have 3-5 parameters that need to be tuned
to reproduce some observed clustering signal, which may decrease the
constraining power of galaxy surveys by introducing degeneracies in the theory
parameter space.

From an observational perspective, an operative definition of bias can be
given in terms of the power spectrum of matter $P_{\rm m}$ and the one of the
biased tracer $P_{\rm t}$ as

$P_{\rm t}(k,\tau)\equiv b_{\rm t}^{2}(k,\tau)P_{\rm m}(k,\tau)\,,$ (1.90)

where the tracer’s bias $b_{\rm t}(k,\tau)$ is a function of both scale and
time.

####  1.3.5 Redshift space distortions

Galaxy surveys provide information on galaxy positions in terms of angular
coordinates and redshift. The redshift we observe can be interpreted as the
result of two effects, the cosmological redshift due to the expansion of the
universe and the Doppler shift due to the peculiar velocity of galaxies in the
reference frame of the expanding universe:

$1+z_{\rm obs}=(1+z_{\rm cos})\cdot\left(1+v_{||}(\mathbf{r})\right)\,,$
(1.91)

where $v_{||}=\mathbf{v}\cdot\hat{r}$ and $\mathbf{r}$ is the vector
originating from the observer and pointing towards the source666Eq. 1.91
assumes $v_{||}\ll 1$ and neglects relativistic effects beyond the first
order. It also assumes that the observer is at rest in the reference frame of
the expanding universe.. In the FLRW metric, the cosmological redshift is due
to the expansion of the universe from the time of emission to that of the
observation

$1+z_{\rm cos}=\frac{a_{\rm obs}}{a_{\rm emit}}\,.$ (1.92)

This can be used to obtain the instantaneous distance from the observer to the
source at the time of the observation

$r=a_{0}\int_{0}^{\chi_{s}}d\chi=a_{0}\int_{t_{s}}^{t_{0}}\frac{dt}{a}=\int_{0}^{z_{\rm
cos}}\frac{dz^{\prime}}{H\left(z^{\prime}\right)}\,,$ (1.93)

where $\chi_{s}$ is the comoving distance of the source. While eq. (1.93)
provides a way to map cosmological redshift to instantaneous positions, the
cosmological redshift of the sources is not directly observed and the observed
redshift is affected also by the peculiar velocities of the sources, as we
discussed above. Using the observed redshift $z_{\rm obs}$ it is possible to
infer the redshift space position

$s\equiv\int_{0}^{z_{\rm
obs}}\frac{dz^{\prime}}{H\left(z^{\prime}\right)}\simeq
r+\frac{v_{||}(\mathbf{r})}{\mathcal{H}\left(z_{\rm cos}\right)}\,.$ (1.94)

On the one hand, this does not allow us to know the exact positions of the
sources, on the other hand, it makes cosmological probes of clustering
sensitive also to the matter velocity field.

In fact, the Doppler shift produces well-studied distortions in the clustering
of galaxies known as Redshift-Space Distortions (RSD). To understand RSD it is
convenient to split the analysis into large (linear) scales and small (non-
linear) scales.

On large scales, the coherent infall of matter towards the overdensity
produces a squashing of the structures along the line-of-sight that is well
described by the Kaiser model in the distant-observer approximation
Kaiser:1987qv. In this model the velocity field is due to the gradient of the
gravitational potential and vorticity can be neglected. This is valid on very
large scales where linear theory applies and the galaxy number density
contrast $\delta_{\rm g}$ is related to the matter density contrast
$\delta_{\rm m}$ by the linear bias relation

$\delta_{\mathrm{g}}(\mathbf{r},\tau)=b(\tau)\delta_{\mathrm{m}}(\mathbf{r},\tau)\,.$
(1.95)

The conservation of galaxy counts under the redshift space remapping
determined by eq. (1.94) relates the redshift-space galaxy number-density
contrast $\delta_{{\rm g},\mathbf{s}}$ to the matter density contrast
$\delta_{\rm m}$ through

$\delta_{\mathrm{g},\mathbf{s}}(\mathbf{r})=b\,\delta_{\mathrm{m}}(\mathbf{r})-\frac{\partial}{\partial\mathbf{r}}\left[\frac{v_{||}(\mathbf{r})}{\mathcal{H}}\right]\,,$
(1.96)

which, in Fourier space, using the Euler equation (1.59) and incorporating
time dependency into the linear growth factor, yields

$\delta_{\mathrm{g},\mathbf{s}}(\mathbf{k})=\left[b+f\mu_{k}^{2}\right]\delta_{\mathrm{m}}(\mathbf{k})\,,$
(1.97)

where $\mu_{\mathbf{k}}$ is the cosine of the angle between the versor
$\hat{k}$ and the line-of-sight, and $f$ is the linear growth rate already
defined in eq. (1.63). Eq. (1.97) shows that the clustering of galaxies
depends on the cosine with the line-of-sight $\mu_{k}$ making redshift-space
galaxies more clustered than real-space galaxies along the line-of-sight. This
clearly affects also the galaxy power spectrum in redshift space which can be
expressed in terms of the matter power spectrum as

$P_{{\rm g},\mathbf{s}}(k,\mu_{k})=\left(b+f\mu_{k}^{2}\right)^{2}P_{\rm
m}(k)\,.$ (1.98)

By expanding the redshift space galaxy power spectrum in its multipoles
through the projection on the Legendre polynomials,

$P_{\ell}^{g,\mathbf{s}}(k)=\frac{2\ell+1}{2}\int_{-1}^{1}P^{g,\mathbf{s}}(k,\mu_{k})\mathcal{L}_{\ell}(\mu_{k})d\mu_{k}\,,$
(1.99)

it is possible to show that, in the Kaiser model, only monopole, quadrupole
and hexadecapole moments are sourced by redshift space distortions, with the
three multipoles given by:

$\displaystyle P_{0}^{g,\mathbf{s}}(k)$
$\displaystyle=\left(b^{2}+\frac{2}{3}f\,b+\frac{1}{5}f^{2}\right)P_{\text{m}}(k)\,,$
(1.100) $\displaystyle P_{2}^{g,\mathbf{s}}(k)$
$\displaystyle=\left(\frac{4}{3}f\,b+\frac{4}{7}f^{2}\right)P_{\text{m}}(k)\,,$
$\displaystyle P_{4}^{g,\mathbf{s}}(k)$
$\displaystyle=\frac{8}{35}f^{2}P_{\text{m}}(k)\,.$

On small scales, instead, the matter is virialised, hence the motion of
galaxies is incoherent and their velocities are larger. The redshift-space
displacement due to the peculiar velocity is larger than the separations
between galaxies which results in an elongation of structures along the line-
of-sight known as the Fingers-of-God effect Jackson:1971sky. Due to their
different nature, unlike the large-scale squashing of structures, these small-
scale distortions require a much more complicated non-linear treatment and, as
a result, they are more often regarded as a problem to deal with rather than a
cosmological probe.

####  1.3.6 Galaxy surveys

Extracting information from the large scale structure is a non-trivial task
and equally (if not more) difficult is to design and run experiments that are
able to accurately map the LSS. It is possible to classify the galaxy surveys
in two main kinds, the ones devoted to collect imaging data (photometric
surveys) and the ones devoted to collect spectrographic data (spectroscopic
surveys). To greater interest of this work are the latter which, unlike the
former, can provide precise determinations for the redshift of the observed
galaxies (necessary to study RSD for example). To achieve this goal in these
experiments first it is necessary to identify a population of galaxies by
means of imaging data, then, known the angular coordinates of each galaxy, the
light collected by the telescope is conveyed to the spectrographs by means of
optical fibers and finally the redshift is determined by analysing the shift
of the observed spectral lines. Depending on the experimental design several
systematic errors need to be taken into account like fiber collisions
Hahn:2016kiy and selection effects 2012MNRAS.424..564R.

Chronologically, the first spectroscopic galaxy survey with significant
cosmological implications is the 2 degrees Field Galaxy Redshift Survey
(2dFGRS) 2DFGRS:2001zay that, from the year 1997 to the 2002, measured the
spectra of more than 200.000 galaxies in a region of $\sim 2.000\,{\rm
deg}^{2}$ around redshift $z=0.11$ 2dFGRS:2005yhx. The analysis of the large
scale galaxy power spectrum in the 2dFGRS provided the first constraints on
the matter content of the universe from LSS observations 2dFGRS:2001csf;
2dFGRS:2005yhx. Soon after the 2dFGRS, other galaxy surveys started collecting
data. The 6 degrees Field Galaxy Survey (6dFGS) carried out between 2001 and
2009 measured the spectra of $\sim 130.000$ galaxies across $\sim 17.000\,{\rm
deg}^{2}$ (almost half of the sky) becoming the widest spectroscopic survey of
the low redshift universe Jones:2004zy; Beutler:2011hx. Between the years 2006
and 2011, the WiggleZ survey Drinkwater:2009sd probed a $\sim 1\,{\rm
Gpc}^{3}$ volume of the universe by focusing on a relatively small area of the
sky ($\sim 1.000\,{\rm deg}^{2}$) but targeting a higher-redshift galaxy
population. With a single robotic spectrograph, this survey collected precise
redshift measurements of about $240.000$ galaxies emission line galaxies in
the redshift range $0.2\lesssim z\lesssim 1$ Blake:2012pj. Across the first 20
years of the 21 century, the Sloan Digital Sky Survey SDSS:2000hjo;
eBOSS:2020yzd measured the spectroscopic redshift of more than 1 million
galaxies in an area of around $\sim 14.000\,{\rm deg}^{2}$, roughly a quarter
of the sky. Having targeted several galaxy populations from nearby galaxies to
Luminous Red Galaxies, Emission Line Galaxies and quasars (or quasi stellar
objects) this survey covered the redshift range $0.07<z<2.2$. With a volume of
more than $10\,{\rm Gpc}^{3}$, SDSS is the largest spectroscopic galaxy
catalogue to date and has enabled precise inference of cosmological
parameters, in particular when combined with imaging surveys like DES and CMB
measurements like Planck.

Motivated by the amazing achievements of these surveys and the relentless
advance in technology and engineering that our times are witnessing, new
surveys have been designed and are now becoming a reality. Among these new
exciting experiments, often referred to as Stage IV surveys, the following are
worthy of a particular mention:

  * •

the Dark Energy Spectroscopic Instrument (DESI) DESI:2016fyo, a ground based
telescope with robotically-orientable optical fibers fed in an array of 10
spectrographs, started collecting data in April 2021 and after only 7 months
of observations became the largest spectroscopic survey available777Source:
Berkeley lab website https://www.lbl.gov.. The experiment is designed to
collect the spectra of magnitude limited bright galaxies at $z\sim 0.2$, LRG
up to $z<1$, ELG up to $z<1.7$, and quasars covering $14.000\,{\rm deg}^{2}$
of the sky. With its wide field and improved sensitivity, DESI will observe a
number of galaxies and a volume $\sim 10$ times larger than SDSS.

  * •

the Euclid mission Euclid:2021icp of the European Space Agency is a space-
based experiment. The launch of the satellite is expected in July 2023 and the
main survey (the wide survey) is planned to start shortly thereafter. Unlike
the others, the Euclid telescope will perform slit-less spectroscopy which may
lead to different systematic errors. The main target of the wide survey will
be the H$\alpha$ emitting galaxies across an area of $15.000\,{\rm deg}^{2}$
in the redshift range $0.9<z<1.8$. Given the estimated number density of these
galaxies Pozzetti:2016cch, Euclid is expected to measure the spectroscopic
redshift of around 30 million galaxies Euclid:2019clj.

  * •

The Nancy Grace Roman space telescope Akeson:2019biv, also known as WFIRST, is
scheduled to launch by May 2027888From the press release: NASA Confirms Roman
Mission’s Flight Design in Milestone Review.. During the first five year of
the mission, WFIRST will devote a significant part of its observing time to
map a deep but relatively narrow region of the sky producing the so-called
High Latitude Survey (HLS). The telescope will measure the spectroscopic
redshifts of galaxies emitting the H$\alpha$ (10 million spectra at z=1-2) and
O-III lines (2 million spectra at z=2-3) in a area of $2.000\,{\rm deg}^{2}$
Wang:2021oec. Similarly to the Euclid mission, WFIRST will perform slit-less
spectroscopy.

As the experiments progress, the theoretical modelling must keep the pace to
maximise the scientific return of Stage IV surveys. Since structure formation
is an intrinsically non-linear process, numerical simulations will play a key
role in future cosmological inference with the LSS.

###  1.4 N-body simulations

Cosmological simulations are based on the N-body framework, where a large
number of particles (typically $\sim 10^{8}-10^{11}$) is used to discretise
the distribution of mass in a comoving portion of the universe referred to as
the simulation box, normally a cubic box with periodic boundary conditions.
They rely on the Newtonian approximation and use comoving coordinates to
incorporate the effects of background expansion. It is also possible to
include large-scales relativistic effects in Newtonian N-body simulations
thanks to the N-body gauge formalism Fidler:2015npa; Tram:2018znz;
Dakin:2019vnj; Brando:2020ouk; Brando:2021jga. The N-body particles are
assumed to be collision-less and to interact only through gravity forces.
Their equations of motions are given by the Hamilton equations

$\displaystyle{\mathbf{x}}^{\prime}_{n}=\frac{\mathbf{p}_{n}}{a(\tau)}$
(1.101)
$\displaystyle{\mathbf{p}}^{\prime}_{n}=a(\tau)\left(-\frac{\partial\phi}{\partial\mathbf{x}}\right)_{n}$
(1.102)

where $\mathbf{x}_{n}$ is the comoving position of the $n^{\rm th}$ particle
and $\mathbf{p}_{n}$ is its conjugate momentum. The gravitational potential
$\phi$ is computed from the density contrast using the Poisson equation
(1.53). The scale factors in the Hamilton equations take into account the
weakening of gravitational interactions due to background expansion. The
number of particles $N_{\rm part}$ and the size of the box $L_{\rm box}$ are
crucial in determining the mass resolution of the simulation. For particles of
equal mass in a cubic box, the mass of each particle is given by

$M_{\rm part}=\bar{\rho}\frac{L_{\mathrm{box}}^{3}}{N_{\mathrm{part}}}\,,$
(1.103)

where $\bar{\rho}$ is the average density of the universe for the given
cosmology.

The Initial Conditions (IC) are normally set using (first or second order)
Lagrangian Perturbation Theory (LPT) in the matter domination at a stage where
nonlinearities are negligible for the scales and accuracy of interest in the
particular case. It has been shown that redshift of $z\gtrsim 100$ should be
used to set the IC with the Zeldovich approximation (i.e., the first order
lagrangian perturbation theory) to achieve percent level accuracy up to $k\sim
1\,h\,{\rm Mpc}^{-1}$ Schneider:2015yka.

The particles’ positions and velocities are evolved from the IC to the desired
redshift using N-body techniques that solve the equations of motion (1.101)
for the N-body particles together with the Poisson equation (1.53) at each
time-step. The solution of the Poisson equation is the most computationally
expensive part in a $\Lambda$CDM simulation, therefore many efforts have been
made to adopt the most efficient techniques. Some of the most relevant
techniques used to compute the forces are:

  * •

Direct-summation method Ewald1921: The equation of motion are integrated once
analytically. The forces are computed by adding the contributions from each
pair of particles, making the computational complexity of this algorithm scale
with $N_{\rm part}^{2}$.

  * •

Particle-Mesh (PM) method: A regular mesh is constructed and the particles are
assigned to the mesh intersections following specific interpolation methods
like Cloud-in-Cell (CIC) or Triangular-Shaped-Clouds (TSC) Sefusatti:2015aex.
This has the advantage that the FFT algorithm can be applied to solve the
Poisson equation in Fourier space. The size of the mesh cells determines the
force resolution $\ell_{\rm F}=L_{\rm box}/N_{\rm mesh}^{1/3}$, where $N_{\rm
mesh}$ is the total number of knots in the mesh. The cosmological code GLAM
Klypin:2017iwu uses the PM method.

  * •

Tree method 1986Natur.324..446B: The particles are assigned to a hierarchical
tree to speed up the force computation. The forces are computed with direct
summation between particle pairs on small scales and between particles and
cells on larger scales. PKDGRAV3 pkdgrav implements the tree method.

  * •

Adaptive Mesh Refinement (AMR) method: Based on the PM method, the AMR method
achieves a better force resolution while mitigating the loss of efficiency by
sampling high-density regions with increasingly finer meshes, until a maximum
refinement level is reached. RAMSES Teyssier:2001cp is an example of AMR
N-body code.

  * •

Tree-PM method: The forces are calculated with the tree method on small scales
and with the PM method on large scales. The codes GADGET3 gadget2; angulo2012
and AREPO Arepo) adopt this technique to carry out the force computation.

The codes PKDGRAV3, RAMSES and GADGET3 have been compared in Schneider:2015yka
for the power spectrum predictions where they have been found to be in $1\%$
agreement up to $k\sim 1\,h\,{\rm Mpc}^{-1}$ with larger discrepancy arising
at smaller scales (in particular at higher redshift). More recently, another
comparison has been carried out using different simulations parameters
Angulo:2020vky finding below percent agreement in the matter power spectrum
between these codes up to $k\sim 5\,h\,{\rm Mpc}^{-1}$ at redshift $z=0$. This
is shown in figure 1.2, depicting a comparison of the matter power spectra at
redshift $z=0$ between different N-body codes. With the exception of
“Gadget3”, relative to the GADGET3 simulation performed in Schneider:2015yka,
all the other results agree are in sub-percent agreement up to $k\sim
5\,h\,{\rm Mpc}^{-1}$. The codes involved in this more recent comparison are
the three also used in Schneider:2015yka (i.e., GADGET3, PKDGRAV3 and RAMSES)
with the addition of ABACUS, a full N-body code based on direct-summation for
near-field forces and multipole expansion for far-field forces.

Fig. 1.2:  Comparison of matter power spectra at redshift $z=0$ obtained with
the N-body codes listed in the legend. Adapted from Angulo:2020vky.

####  1.4.1 N-body simulations in modified gravity

When MG is taken into account, the Poisson equation is modified accordingly
(see eq. (1.13)) and the Klein-Gordon equation for the scalar field also need
to be solved at each time step. The latter is a highly non-linear partial
differential equation that needs to be solved with multi-grid techniques,
where the equation is discretised on a grid (possibly using AMR) and solved
iteratively with the Gauss-Seidel scheme Llinares:2018maz. This makes
cosmological simulations in MG much more computationally expensive than their
$\Lambda$CDM counterparts, being typically $\sim 10$ times slower
Winther:2014cia.

Some of the N-body codes in MG that solve the full differential equation for
the scalar field have been compared in Winther:2015wla. The codes included in
the comparison are:

  * •

ECOSMOG Li:2011vk; Li:2013nua, based on the AMR code RAMSES

  * •

MG-GADGET Puchwein:2013lza, based on the tree-PM code GADGET3, with AMR for
the MG solver

  * •

DGPM Schmidt:2009sg; Oyaizu:2008sr, fixed grid PM code

  * •

ISIS Llinares:2013jza, also based on the AMR code RAMSES

  * •

ISIS-NONSTATIC Llinares:2013qbh, version of ISIS that goes beyond the quasi-
static approximation

ECOSMOG, MG-GADGET and ISIS are used to produce $f(R)$ simulations. DGPM and
ECOSMOG are used to produce simulations in DGP theories999The comparison in
DGP theories is performed using a fixed PM grid in ECOSMOG to match the force
resolution of DGPM which does not incorporate an AMR solver.. Finally, ISIS
and ISIS-NONSTATIC are used to produce simulations in a Symmetron model. The
codes are compared at the level of the matter power spectrum, velocity
divergence, halo mass function and halo profiles. Since differences in these
statistics between the various codes are already present in GR, the comparison
focuses on the enhancement of the statistics in MG compared to GR (the MG
boost factor). For what concerns the MG boost factor of the power spectrum at
redshift $z=0$, all codes are found to be in $1\%$ agreement up to $k\sim
7\,h\,{\rm Mpc}^{-1}$ in $f(R)$ and DGP theories. Also, the agreement in the
velocity divergence power spectrum is found to be excellent, with better than
$1\%$ agreement for $k\lesssim 3h\mathrm{Mpc}^{-1}$ in $f(R)$ theory and
better than $\lesssim 2\%$ at all scales considered for DGP theories. The halo
statistics are found to be in good agreement too in both $f(R)$ and DGP
theories, but sample variance limits the validity of the comparison.
Interestingly, ISIS and ISIS-NONSTATIC are found to give almost identical
results for all the statistics considered in the symmetron model, suggesting
that it is fine to neglect the time derivatives in the Klein-Gordon equation
with the quasi-static approximation.

More recently, other simulation codes have been extended to MG. Of interest to
this work are the MG extension of AREPO and GLAM, named MG-AREPO Arepo_fR;
Arepo_nDGP and MG-GLAM Hernandez-Aguayo:2021kuh; Ruan:2021wup, that we will
encounter again in chapter 4.

An alternative way to incorporate screening in N-body simulations in MG is
thanks to the screening approximation proposed in Winther:2014cia and that we
have discussed in the case of chameleon and Vainshtein mechanisms in section
1.2.1. This approximation, which can be obtained by linearising the Klein-
Gordon equation, allows reducing the slow-down of N-body simulations due to MG
to as little as a factor $2$-$3$ while still capturing the suppression of the
fifth force on small scales.

##  Chapter 2 Fast production of galaxy mock catalogues in modified gravity

> The content of this chapter is based on the publication Fiorini:2021dzs. The
> particles data from ELEPHANT simulations used in this chapter were provided
> by Baojiu Li.

It is crucial to make accurate theoretical predictions of the properties of
the LSS on non-linear scales, both in terms of comparison to measurements and
to construct realistic covariance matrices, to successfully constrain the
gravity model (see e.g. Alam:2020jdv). This can be achieved by means of
cosmological simulations, but in excess of ${\cal O}(10^{3})$ realisations
must be produced to match the volume of Stage IV surveys and compute an
accurate estimate of the covariance matrices. This poses a serious challenge
for full N-body simulations in MG models due to their high computational cost.
Such models usually have screening mechanisms that hide modified gravity
effects on small scales. To describe the screening mechanism, an additional
non-linear equation needs to be solved in N-body simulations, which
significantly slows down the MG simulations, as discussed in section 1.2.

An alternative to full N-body simulations is to exploit approximate methods
(see Monaco16 for a comprehensive review and Chuang:2014toa for a comparison
project in GR), which lower the computational cost at the expenses of accuracy
on non-linear scales. Amongst these, the COLA method Tassev:2013pn;
Koda:2015mca; Izard:2015dja; Howlett:2015hfa with its extension to MG
Valogiannis:2016ane; Winther:2017jof; Wright:2017dkw offers an interesting
compromise between speed-up and accuracy without introducing any additional
free parameter and is therefore ideal to access the MG information on (mildly)
non-linear scales without losing predictability (see Moretti:2019bob for an
alternative approach based on pinocchio111pinocchio uses analytical LPT
solutions to evolve the DM distribution and the orbit-crossing collapse model
to produce halo catalogues Monaco:2001jg.).

Having a bridge between observed galaxies and simulated dark matter
distribution makes it possible to extract the LSS information encoded in the
datasets from galaxy surveys. As we have seen in section 1.3.4, the HOD
prescription allows us to do so starting from the DM halos identified in the
density field by a halo-finder algorithm Knebe:2011rx and populating them with
galaxies with a probability distribution conditioned on some halo properties.
In particular, the HOD model proposed in Zheng:2007zg relies on the halo mass
and the density: using the analytical NFW model (introduced in sub-section
2.2) for the halo profile, it is possible to apply this formalism to COLA
simulations that do not resolve the internal halo properties Koda:2015mca. Due
to the non-trivial dynamics of the screening mechanisms, the internal halo
structure in MG theories can differ from the one in GR and this needs to be
considered to produce realistic galaxy mocks Mitchell:2018qrg;
Mitchell:2019qke.

In this context, we investigate the feasibility of producing mock galaxy
catalogues from COLA simulations in MG producing a simulation suite with mg-
picola Winther:2015wla; Wright:2017dkw and using a suite of full N-body
simulations performed by the ecosmog code Li:2011vk to validate the COLA
results and to estimate the accuracy in reproducing galaxy clustering
statistics.

This chapter is organised as follows. After introducing the MG theories of
interest to this work in Section 1.2, we discuss the simulation techniques and
suites in Section 2.1. In Section 2.2 we investigate the production of halo
catalogues. We then apply the HOD formalism in Section 2.3 to create mock
galaxy catalogues and study the multipole moments of the galaxy spectrum in
redshift space.

###  2.1 Simulations

To produce mock catalogues for the modified gravity models previously
described, we start by running dark matter simulations of the large-scale
structure. In particular, we explore the following modified gravity models:
the normal branch DGP model with $H_{0}r_{c}=1$ (N1) and Hu-Sawicki $f(R)$
model with $|f_{R0}|=10^{-5}$ (F5), plus the vanilla GR model.

Cosmological simulations can be accurate up to non-linear scales, but they are
computationally very expensive, as forces between particles need to be
calculated for thousands of time steps typically. The cost is even higher for
MG models, which in addition also need to solve the non-linear equation for
the fifth force. In recent years, the so-called approximate methods have
become a popular alternative to full $N$-body simulations, and they are aimed
at speeding up the creation of a density field at the expense of not resolving
accurately sub-halo scales (see Chuang:2014toa; Lippich_2018; Blot_2019;
Colavincenzo_2018 for studies on their accuracy, and Monaco16 for a review).
The simulations in this work make use of the COmoving Lagrangian Acceleration
(COLA) method Tassev:2013pn, in particular, we use mg-picola Winther:2017jof,
that includes gravity models other than GR, and is an ideal tool to
efficiently run simulations for MG.

####  2.1.1 COLA method

The COLA method Tassev:2013pn uses the Particle-Mesh (PM) algorithm (see sec
1.4) to evolve the displacement of particles with respect to their second-
order Lagrangian Perturbation Theory (2LPT) positions which are obtained
analytically (see sec 1.3.2). In practice, the PM technique in COLA solves the
following equation of motion for the particles trajectories,

$\partial_{t}^{2}\mathbf{x}_{\mathrm{res}}=-\mathbf{\nabla}\Phi-\partial_{t}^{2}\mathbf{x}_{\mathrm{2LPT}}\,,$
(2.1)

where $\mathbf{x}_{\mathrm{2LPT}}=q+\Psi^{(1)}+\Psi^{(2)}$ is the 2LPT
solution that is computed beforehand. The second time derivative of
$\mathbf{x}_{\mathrm{2LPT}}$ acts as a fictitious force in the (non-inertial)
COLA reference frame. Thanks to this, a small number of time-steps (typically
of $\mathcal{O}(10)$) are enough to accurately recover the DM density field up
to mildly non-linear scales, as well as the halo positions and masses. This
results in a speed-up of a factor $100-1000$ with respect to full N-body
simulations, at the expense of not resolving internal halo properties because
of the low force and time resolution that are used for the PM technique. For
the production of accurate mock halo catalogues, Izard:2015dja proposed an
optimal configuration of the parameters in COLA that control the trade-off
between accuracy and computational cost, such as the grid resolution used to
compute forces and the distribution of time steps, which we adopt for this
work.

####  2.1.2 MG extension to COLA

When MG is incorporated in cosmological simulation, the complexity increases
because the motion of particles is governed by both the Newtonian force and
the fifth force. This can lead to a slow-down of a factor of $\mathcal{O}(10)$
in theories where computing the fifth force requires the solution of non-
linear equations. One way to avoid this significant slow-down is to use the
screening approximations introduced in section 1.2 where the fifth force
dynamics is described by means of an effective mass of the scalar field and a
coupling determined by either the value of gravitational potential or its
derivatives, depending on the type of screening mechanism. To extend the COLA
method to MG, it is also necessary to reformulate the Lagrangian Perturbation
Theory to take into account the effect of the fifth force on the dynamics
Valogiannis:2016ane; Winther:2017jof; Aviles:2017aor.

mg-picola Winther:2017jof is a publicly available code for cosmological
simulations with the COLA method in modified gravity and it implements the
scale-dependent 2LPT and some screening approximations. Its accuracy has
already been tested against full N-body simulations at the level of the DM
density field. In this work, we study the accuracy of mg-picola also for the
DM velocity field, the halo abundance and halo clustering, and we show how it
can be used in combination with a HOD algorithm to generate mock galaxy
catalogues (with some tweaks relative to applications in GR).

####  2.1.3 Simulation suites

We use a full N-body simulation suite called elephant, which was introduced
and validated in Cautun:2017tkc, to benchmark the results of our runs based on
COLA. The elephant suite was produced using ecosmog Li:2011vk, a MG extension
of the adaptive mesh refinement code ramses Teyssier:2001cp that solves the
exact equations for the fifth force; these simulations were recently used to
create mock galaxy catalogues in Hernandez-Aguayo:2018yrp; Hernandez-
Aguayo:2018oxg; Alam:2020jdv222see Barreira:2016ovx and Devi:2019swk for other
attempts to create mock galaxy catalogues in MG theories. Table 2.2 describes
the parameters of the suite such as the box size, the mass resolution and the
gravity models implemented, and the cosmological parameters employed are:

$\begin{array}[]{ccc}\Omega_{m,0}=0.281\,,&\Omega_{\Lambda,0}=0.719\,,&\Omega_{b,0}=0.046\,,\\\
n_{s}=0.971\,,&\sigma_{8}=0.842\,,&h=0.697\,.\end{array}$ (2.2)

The initial conditions were generated using the Zel’dovich approximation at
$z=49$. For a given realisation, the same initial seed was used for GR, F5 and
N1333Due to an anomaly with one snapshot, we use only 4 realisations for N1 in
this work..

We develop a new suite of simulations with mg-picola that matches the
cosmology and parameters (such as the mass resolution, the box size and the
number of realisations) of the fiducial elephant set. We call this new COLA
suite PIpeline TEsting Run (piter). We stop the simulations at redshift
$z=0.5057$ after 30 timesteps to match the redshift of the elephant snapshot
closest to redshift $z=0.5$. As discussed in the previous section, with the
screening approximation it is possible to tune the behaviour of the screening
mechanism: in Eq. (1.40) we use a screening efficiency $f_{s}=4.0$ for F5 and
in Eq. (1.44) we use a Gaussian smoothing with scale 1 $\,h^{-1}\,{\rm Mpc}$
for N1. We have found these values to provide the best agreement for the dark
matter power spectrum with N-body simulations amongst the values tested. A
summary of the piter simulation details is given in Table 2.2.

Table 2.1: elephant simulations Models | GR, F5, N1
---|---
Realisations | 5
Box size |  1024 $\,h^{-1}\,{\rm Mpc}$
$N_{\mathrm{part}}$ | $1024^{3}$
$M_{part}$ | $7.7985\cdot 10^{10}\,h^{-1}\,M_{\odot}$
Domain grid | $1024^{3}$
Refinement criterion | $8$
Initial conditions |  Zel’dovich, $z=49$

Table 2.2: piter simulations Models | GR, F5, N1
---|---
Realisations | 5
Box size |  1024 $\,h^{-1}\,{\rm Mpc}$
$N_{\mathrm{part}}$ | $1024^{3}$
$M_{part}$ | $7.7985\cdot 10^{10}\,h^{-1}\,M_{\odot}$
Force grid | $3072^{3}$
Timesteps | $30$
Initial conditions |  2LPT, $z=49$

Initial conditions are generated using the 2LPT at $z=49$. As in elephant, for
a given realisation, the same initial seed is used for GR, F5 and N1. We also
ran a GR COLA simulation using the same initial condition as one of the
realisations in elephant and checked that our conclusions were not affected by
the cosmic variance. For simplicity of notation, in the following we will just
use N-body and COLA to refer to elephant and piter respectively.

We use the public codes PowerI4444https://github.com/sefusatti/PowerI4 and
DTFE Cautun:2011gf to compute and compare the density and velocity-divergence
power spectra in the elephant and piter simulations. In Figures 2.1 and 2.2 we
compare the clustering signals in GR (top left panel), the ratio between COLA
and N-body (top right panel) and the enhancement in MG with respect to GR,
which we refer as the boost-factor (bottom panels, F5 on the left and N1 on
the right). The measurements correspond to the mean over 5 realisations and
the shaded regions display the standard deviation (just to have a rough
estimate of the errors to guide the eye). The errors are smaller in the boost-
factor plots because the noise cancels out: the MG and the GR simulations were
started from the same initial conditions and therefore have a very similar
cosmic variance at late redshifts that vanishes in ratios. We find better than
$3\%$ agreement between COLA and N-body for the density power spectra up to
modes of $\sim 1\,h\,{\rm Mpc}^{-1}$ in GR while the boost-factors in F5 and
N1 show agreement within the variance up to $k\sim 1\,h\,{\rm Mpc}^{-1}$. This
shows that COLA simulations capture the MG effects very accurately even with
the use of approximations for screening. The agreement in the velocity
divergence power spectra is within $3\%$ up to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$
in GR and for the boost factor in F5 while the boost factor in N1 agrees
within $2\%$ up to $k\sim 2\,h\,{\rm Mpc}^{-1}$. In GR, after losing power at
$k\sim 1\,h\,{\rm Mpc}^{-1}$, the velocity divergence power spectrum in COLA
shows an enhancement compared with N-body at very high $k$. This is because
COLA has less random motion at small scales that washes out the signal in
N-body more than in COLA. The better agreement in the density power spectra is
expected because the velocity divergence power spectrum is more affected by
non-linear physics.

From the boost-factors, we can already appreciate the intrinsic difference
between F5 and N1: while F5 has a scale-dependent enhancement, N1 is almost
scale-invariant across most scales. This behavior can be explained by
comparing the respective Poisson equations: $G_{\rm eff}$ is a function of $k$
in Eq. (1.37) for $f(R)$ gravity while it is scale-independent in Eq. (1.43)
for DGP.

(a) GR

(b) GR ratio

(c) F5 boost factor

(d) N1 boost factor

Fig. 2.1:  Dark Matter power spectrum comparison between COLA (orange dashed
lines) and N-body (blue solid lines) obtained by taking the average over 5
realisations. The shaded regions represent the standard deviation over the 5
realisations. The power spectrum in GR (top panels) and the boost-factors in
F5 (bottom left) and N1 (bottom right) show agreement within the variance up
to $k\sim 1\,h\,{\rm Mpc}^{-1}$.

(a) GR

(b) GR ratio

(c) F5 boost factor

(d) N1 boost factor

Fig. 2.2:  Dark Matter velocity divergence power spectrum comparison between
COLA (orange dashed lines) and N-body (blue solid lines) obtained by taking
the average over 5 realisations. The shaded regions represent the standard
deviation over the 5 realisations. The power spectrum in GR (top panels) is in
agreement within $3\%$ up to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$ where COLA starts
to lose power. The boost-factor in F5 (bottom left) and N1 (bottom right) show
better than $5\%$ accuracy even deeper in the non-linear regime, up to $k\sim
2\,h\,{\rm Mpc}^{-1}$.

###  2.2 Halos

Halos can be found by means of specific algorithms often referred to as halo-
finders. Amongst the several halo-finders proposed in literature we focus on
two popular codes: rockstar Behroozi13 and the Friends-of-Friends (FoF) finder
Davis85 (as implemented in nbodykit555https://nbodykit.readthedocs.io).
rockstar uses a 6D phase-space FoF finder and is capable of measuring many
halo properties like the spherical over-density mass, the concentration
parameter and the velocity dispersion just to name a few. The FoF finder
instead is purely based on particles positions.

####  2.2.1 rockstar halos

We first produce rockstar halo catalogues in COLA and N-body and compare the
cumulative halo mass function (hmf) and the halo power spectrum. We use the
setting to remove unbound particles and keep only parent halos666We found
that, if we do not remove unbound particles, the agreement between N-body and
COLA is substantially worse for the halo power spectrum. The hmf is simply
defined as the number density of halos with mass greater than $M$. Before
measuring the halo power spectrum, we draw a halo sample defined by an
abundance matching procedure, in which the target density is given by a cut at
$10^{13}\,h^{-1}\,M_{\odot}$ in the N-body catalogues for GR. Figures 2.3 and
2.4 show the comparisons of the hmf and the halo power spectrum for COLA and
N-body rockstar catalogues respectively, both in GR. The hmf is $25\%$ off for
masses greater than $10^{13}\,h^{-1}\,M_{\odot}$, and the halo power spectrum
ratio is well predicted in the large-scale limit but shows a scale-dependent
deviation, up to $10\%$ at $k\sim 0.4\,h\,{\rm Mpc}^{-1}$. These differences
are likely due to the lower force resolution and accuracy in recovering the
velocity field in COLA simulations.

(a) GR

(b) GR ratio

Fig. 2.3: rockstar halo mass function comparison in GR. The hmf in COLA and
N-body (left panel) and their ratio (right panel) show $\sim 25\%$ discrepancy
for masses greater than $10^{13}\,h^{-1}\,M_{\odot}$.

(a) GR

(b) GR ratio

Fig. 2.4: rockstar power spectrum comparison between the power spectra of COLA
and N-body in GR. The discrepancy, which is already clear in the left panel,
is shown to reach the $10\%$ with a scale-dependent behaviour in the right
panel.

Given this result, in the rest of the chapter, we will use only the FoF finder
to generate halo catalogues from COLA and N-body simulations. However, we
still use rockstar in the N-body simulations to quantify the effects of MG on
the concentration and velocity dispersion of halos in F5 simulations (see
Section 2.2.3), which we then incorporate in the prescriptions that we use to
populate the COLA simulations with galaxies.

####  2.2.2 Friends-of-Friends halos

The FoF algorithm links together particles that are separated by a distance
smaller than a certain linking length $b$, expressed in units of the mean
inter-particle distance. We adopt the standard value $b=0.2$. This value,
which is commonly used in N-body simulations, has been shown to also be valid
in COLA simulations (see Koda:2015mca; Howlett:2014opa; Izard:2015dja),
although, in some other approximate methods that use very low-resolution
density fields, $b$ has been tweaked, see Manera:2012sc. At this point, we are
interested in converting FoF masses, $M_{\rm FoF}$, to spherical over-density
(SO) masses, $M_{\rm 200c}$, which we adopt in the rest of the chapter when
expressing halo masses. $M_{\rm 200c}$ denotes the mass enclosed in a
spherical region with density 200 times the critical density. We use the fit
in Lukic:2008ds to convert mass definitions with an accuracy of $5\%$ for most
halos in the mass range $10^{12.5}-10^{15.5}\,h^{-1}\,M_{\odot}$:

$\frac{M_{\rm
FoF}}{M_{200\mathrm{c}}}=\frac{a_{1}}{c_{200}^{2}}+\frac{a_{2}}{c_{200}}+a_{3}\,,$
(2.3)

where the coefficients $a_{1}$, $a_{2}$ and $a_{3}$ depend on the number of
particles (see their Table 1), and $c_{200}$ is the concentration parameter,
which we compute using the empirical formula from Dutton:2014xda:

$\displaystyle\log_{10}c_{200}$ $\displaystyle=a+b\log_{10}\left(\frac{M_{\rm
200c}}{10^{12}\,h^{-1}\,M_{\odot}}\right),$ (2.4) $\displaystyle a$
$\displaystyle=0.520+(0.905-0.520)\exp\left(-0.617z^{1.21}\right),$
$\displaystyle b$ $\displaystyle=-0.101+0.026z.$

The empirical relations in eq. (2.3) and eq. (2.4) have been calibrated in
$\Lambda$CDM simulations but given the small amplitude of MG effects on the
halo concentration that we discuss in section 2.2.3, we argue that the above
relations are valid within our target accuracy also in the MG theories we
consider in this chapter.

The left panel in Figure 2.5 shows the conversion factor of Eq. (2.3) (for a
set of values corresponding to SO masses of 100, 600, 1000, 3000, 6000 and
10000 particles). There is very little dependence on halo mass, with all the
values within $2\%$ of 0.75. Hence we use the constant value 0.75 to convert
the FoF masses to SO masses. To validate this mass conversion, we compare the
mass function for the FoF halo catalogues $n_{\rm FoF}$ with rescaled
$M_{200\mathrm{c}}=0.75M_{\rm FoF}$ to the mass function as measured by
rockstar $n_{\rm Rock}$ in terms of $M_{\rm 200c}$. The right panel of Figure
2.5 shows this comparison for N-body GR halo catalogues, and the agreement of
these two mass functions is within $3\%$ for masses grater than
$10^{13}\,h^{-1}\,M_{\odot}$. We thus produce halo catalogues for both the
N-body and COLA simulations using the FoF method and assign the over-density
mass $M_{200\mathrm{c}}=0.75M_{\rm FoF}$.

(a) Mass conversion factor

(b) rockstar comparison

Fig. 2.5: Mass conversion between FoF and $M_{200c}$ spherical over-density
masses. The conversion factor is evaluated for 6 mass values (left panel),
corresponding to the SO mass of halos composed of 100, 600, 1000, 3000, 6000
and 10000 DM particles. All the six points lie within $2\%$ from the value