Without sufficiently accurate analytic models of matter clustering, simulations can instead act as a model of nonlinear dynamics.
N-body simulations provide Monte Carlo realizations of Newtonian dynamics in the fully-nonlinear regime, producing a nonlinear matter density field for a given cosmological model.
As with PT and the halo model, N-body simulations are limited by a fundamental assumption - namely that evolution of the matter distribution is fully described by collisionless cold dark matter obeying Newtonian gravity.
However, this assumption does not limit the scales accessible to the model or the types of nonlinear structures it can produce, which are constrained only by numerical resolution.
Recent advances in computing have led to the rise of large-volume, high-resolution simulations (e.g. [1]}, [2]}, [3]}, [4]}, [5]}), albeit with questions of convergence at the percent-level [6]}.
With these have come approximate methods of simulation that aim to obtain comparable solutions with much less computation time [7]}, [8]}, [9]}.
As redshift surveys push to larger volumes and higher number density, the intractability of running many sufficiently resolved simulations at multi-Gpc volumes (necessary for capturing large modes and estimating covariance) to produce two-point statistics has motivated various fitting functions and interpolations of two-point statistics produced by high resolution simulations (e.g. [10]}, [11]}, [12]}, [13]}, [14]}).
At the smallest scales probed by observations, baryonic effects on the matter distribution are also a concern [15]}, [16]}, and to properly simulate their impact on large-scale structure requires a full understanding of galaxy formation and hydrodynamic simulations that include feedback [17]}, [18]}.
There have been some recent efforts to correct for baryons and mitigate this issue by modifying the output of dark-matter-only simulations [19]}, [20]}, [21]}.
