This is the first in a series of papers devoted to fully general-relativistic $N$-body simulations applied to late-time cosmology. The purpose of this paper is to present the combination of a numerical relativity scheme, discretization method and time-integration algorithm that provides satisfyingly stable evolution. More precisely, we show that it is able to pass a robustness test and to follow scalar linear modes around an expanding homogeneous and isotropic space-time. Most importantly, it is able to evolve typical cosmological initial conditions on comoving scales down to tenths of megaparsecs with controlled constraint and energy-momentum conservation violations all the way down to the regime of strong inhomogeneity.

We present the new relativistic cosmological particle-mesh code $\texttt{fRevolution}$, based on $\texttt{gevolution}$, aimed at simulating non-linear structure formation in $f(R)$ gravity. We introduce the general framework and approximation scheme, and the set of equations used to solve for the full set of gravitational perturbations. We show results for a point mass field and for cosmological simulations in the Hu-Sawicki model, and compare them to those of existing Newtonian codes. A more detailed analysis and discussion of our solutions will be carried out in a following paper.

We present results from the first simulations of networks of Type I Abelian Higgs cosmic strings to include both matter and radiation eras and Cosmic Microwave Background (CMB) constraints. In Type I strings, the string tension is a slowly decreasing function of the ratio of the scalar and gauge mass-squared, $\beta$. We find that the mean string separation shows no dependence on $\beta$, and that the energy-momentum tensor correlators decrease approximately in proportion to the square of the string tension, with additional O(1) correction factors which asymptote to constants below $\beta \lesssim 0.01$. Strings in models with low self-couplings can therefore satisfy current CMB bounds at higher symmetry-breaking scales. This is particularly relevant for models where the gauge symmetry is broken in a supersymmetric flat direction, for which the effective self-coupling can be extremely small. If our results can be extrapolated to $\beta \simeq 10^{-15}$, even strings formed at $10^{16}$ GeV (approximately the grand unification scale in supersymmetric extensions of the Standard Model) can be compatible with CMB constraints.

A persistent theme in the study of dark energy is the question of whether it really exists or not. It is often claimed hat we are mis-calculating the cosmological model by neglecting the effects associated with averaging over large-scale structures. In the Newtonian approximation this is clear: there is no effect. Within the full relativistic picture this remains an important open question however, owing to the complex mathematics involved. We study this issue using particle numerical simulations which account for all relevant relativistic effects without any problems from shell crossing. In this context we show for the first time that the backreaction from structure can differ by many orders of magnitude depending upon the slicing of spacetime one chooses to average over. In the worst case, where smoothing is carried out in synchronous spatial surfaces, the corrections can reach ten percent and more. However, when smoothing on the constant time hypersurface of the Newtonian gauge, backreaction contributions remain 3-5 orders of magnitude smaller.

Cosmological simulations involving the fully covariant gravitational dynamics may prove relevant in understanding relativistic/non-linear features and, therefore, in taking better advantage of the upcoming large scale structure survey data. We propose a new 3+1 integration scheme for General Relativity in the case where the matter sector contains a minimally-coupled perfect fluid field. The original feature is that we completely eliminate the fluid components through the constraint equations, thus remaining with a set of unconstrained evolution equations for the rest of the fields. This procedure does not constrain the lapse function and shift vector, so it holds in arbitrary gauge and also works for arbitrary equation of state. An important advantage of this scheme is that it allows one to define and pass an adaptation of the robustness test to the cosmological context, at least in the case of pressureless perfect fluid matter, which is the relevant one for late-time cosmology.

We investigate cosmic string networks in the Abelian Higgs model using data from a campaign of large-scale numerical simulations on lattices of up to $4096^3$ grid points. We observe scaling or self-similarity of the networks over a wide range of scales, and estimate the asymptotic values of the mean string separation in horizon length units $\dot{\xi}$ and of the mean square string velocity $\bar v^2$ in the continuum and large time limits. The scaling occurs because the strings lose energy into classical radiation of the scalar and gauge fields of the Abelian Higgs model. We quantify the energy loss with a dimensionless radiative efficiency parameter, and show that it does not vary significantly with lattice spacing or string separation. This implies that the radiative energy loss underlying the scaling behaviour is not a lattice artefact, and justifies the extrapolation of measured network properties to large times for computations of cosmological perturbations. We also show that the core growth method, which increases the defect core width with time to extend the dynamic range of simulations, does not introduce significant systematic error. We compare $\dot{\xi}$ and $\bar v^2$ to values measured in simulations using the Nambu-Goto approximation, finding that the latter underestimate the mean string separation by about 25%, and overestimate $\bar v^2$ by about 10%. The scaling of the string separation implies that string loops decay by the emission of massive radiation within a Hubble time in field theory simulations, in contrast to the Nambu-Goto scenario which neglects this energy loss mechanism. String loops surviving for only one Hubble time emit much less gravitational radiation than in the Nambu-Goto scenario, and are consequently subject to much weaker gravitational wave constraints on their tension.

We present cosmic microwave background (CMB) power spectra from recent numerical simulations of cosmic strings in the Abelian Higgs model and compare them to CMB power spectra measured by Planck. We obtain revised constraints on the cosmic string tension parameter $G\mu$. For example, in the $\Lambda$CDM model with the addition of strings and no primordial tensor perturbations, we find $G\mu < 2.0 \times 10^{-7}$ at 95% confidence, about 20% lower than the value obtained from previous simulations, which had 1/64 of the spatial volume. We investigate the source of the difference, showing that the main cause is an improved treatment of the string evolution across the radiation-matter transition. The increased computational volume also makes possible to simulate fully the physical equations of motion, in which the string cores shrink in comoving coordinates. This, and the larger dynamic range, changes the amplitude of the power spectra by only about 10%, demonstrating that field theory simulations of cosmic strings have now reached the required dynamic range for CMB calculations.

We report on the energy-momentum correlators obtained with recent numerical simulations of the Abelian Higgs model, essential for the computation of cosmic microwave background and matter perturbations of cosmic strings. Due to significant improvements both in raw computing power and in our parallel simulation framework, the dynamical range of the simulations has increased four-fold both in space and time, and for the first time we are able to simulate strings with a constant physical width in both the radiation and matter eras. The new simulations improve the accuracy of the measurements of the correlation functions at the horizon scale and confirm the shape around the peak. The normalization is slightly higher in the high wave-number tails, due to a small increase in the string density. We study for the first time the behaviour of the correlators across cosmological transitions, and discover that the correlation functions evolve adiabatically, ie the network adapts quickly to changes in the expansion rate. We propose a new method for constructing source functions for Einstein-Boltzmann integrators, comparing it with two other methods previously used. The new method is more consistent, easier to implement, and significantly more accurate.

We present a new N-body code, gevolution, for the evolution of large scale structure in the Universe. Our code is based on a weak field expansion of General Relativity and calculates all six metric degrees of freedom in Poisson gauge. N-body particles are evolved by solving the geodesic equation which we write in terms of a canonical momentum such that it remains valid also for relativistic particles. We validate the code by considering the Schwarzschild solution and, in the Newtonian limit, by comparing with the Newtonian N-body codes Gadget-2 and RAMSES. We then proceed with a simulation of large scale structure in a Universe with massive neutrinos where we study the gravitational slip induced by the neutrino shear stress. The code can be extended to include different kinds of dark energy or modified gravity models and going beyond the usually adopted quasi-static approximation. Our code is publicly available.

latfield2 is a C++ library designed to simplify writing parallel codes for solving partial differen- tial equations, developed for application to classical field theories in particle physics and cosmology. It is a significant rewrite of the latfield framework, moving from a slab domain decomposition to a rod decomposition, where the last two dimension of the lattice are scattered into a two dimensional process grid. Parallelism is implemented using the Message Passing Interface (MPI) standard, and hidden in the basic objects of grid-based simulations: Lattice, Site and Field. It comes with an integrated parallel fast Fourier transform, and I/O server class permitting computation to continue during the writing of large files to disk. latfield2 has been used for production runs on tens of thousands of processor elements, and is expected to be scalable to hundreds of thousands.

Numerical simulations are a versatile tool providing insight into the complicated process of structure formation in cosmology. This process is mainly governed by gravity, which is the dominant force on large scales. To date, a century after the formulation of general relativity, numerical codes for structure formation still employ Newton's law of gravitation. This approximation relies on the two assumptions that gravitational fields are weak and that they are only sourced by non-relativistic matter. While the former appears well justified on cosmological scales, the latter imposes restrictions on the nature of the "dark" components of the Universe (dark matter and dark energy) which are, however, poorly understood. Here we present the first simulations of cosmic structure formation using equations consistently derived from general relativity. We study in detail the small relativistic effects for a standard {\Lambda}CDM cosmology which cannot be obtained within a purely Newtonian framework. Our particle-mesh N-body code computes all six degrees of freedom of the metric and consistently solves the geodesic equation for particles, taking into account the relativistic potentials and the frame-dragging force. This conceptually clean approach is very general and can be applied to various settings where the Newtonian approximation fails or becomes inaccurate, ranging from simulations of models with dynamical dark energy or warm/hot dark matter to core collapse supernova explosions.

In this work we discuss the possibility of cosmic defects being responsible for the B-mode signal measured by the BICEP2 collaboration. We also allow for the presence of other cosmological sources of B-modes such as inflationary gravitational waves and polarized dust foregrounds, which might contribute to or dominate the signal. On the one hand, we find that defects alone give a poor fit to the data points. On the other, we find that defects help to improve the fit at higher multipoles when they are considered alongside inflationary gravitational waves or polarized dust. Finally, we derive new defect constraints from models combining defects and dust. This proceeding is based on previous works [1,2].

We analyse the possible contribution of topological defects to cosmic microwave anisotropies, both temperature and polarisation. We allow for the presence of both inflationary scalars and tensors, and of polarised dust foregrounds that may contribute to or dominate the B-mode polarisation signal. We confirm and quantify our previous statements that topological defects on their own are a poor fit to the B-mode signal. However, adding topological defects to a models with a tensor component or a dust component improves the fit around $\ell=200$. Fitting simultaneously to both temperature and polarisation data, we find that textures fit almost as well as tensors ($\Delta\chi^2 = 2.0$), while Abelian Higgs strings are ruled out as the sole source of the B-mode signal at low $\ell$. The 95% confidence upper limits on models combining defects and dust are $G\mu < 2.7\times 10^{-7}$ (Abelian Higgs strings), $G\mu < 9.8\times 10^{-7}$ (semilocal strings) and $G\mu < 7.3\times 10^{-7}$ (textures), a small reduction on the $Planck$ bounds. The most economical fit overall is obtained by the standard $\Lambda$CDM model with a polarised dust component.

We show that the B-mode polarization signal detected at low multipoles by BICEP2 cannot be entirely due to topological defects. This would be incompatible with the high-multipole B-mode polarization data and also with existing temperature anisotropy data. Adding cosmic strings to a model with tensors, we find that B-modes on their own provide a comparable limit on the defects to that already coming from Planck satellite temperature data. We note that strings at this limit give a modest improvement to the best-fit of the B-mode data, at a somewhat lower tensor-to-scalar ratio of $r \simeq 0.15$.

We develop a formalism for General Relativistic N-body simulations in the weak field regime, suitable for cosmological applications. The problem is kept tractable by retaining the metric perturbations to first order, the first derivatives to second order and second derivatives to all orders, thus taking into account the most important nonlinear effects of Einstein gravity. It is also expected that any significant "backreaction" should appear at this order. We show that the simulation scheme is feasible in practice by implementing it for a plane-symmetric situation and running two test cases, one with only cold dark matter, and one which also includes a cosmological constant. For these plane-symmetric situations, the deviations from the usual Newtonian N-body simulations remain small and, apart from a non-trivial correction to the background, can be accurately estimated within the Newtonian framework. The correction to the background scale factor, which is a genuine "backreaction" effect, can be robustly obtained with our algorithm. Our numerical approach is also naturally suited for the inclusion of extra relativistic fields and thus for dark energy or modified gravity simulations.