This review is an uptodate account of the use of numerical relativity to study dynamical stronggravity environments in a cosmological context. First we provide a gentle introduction into the use of numerical relativity in solving cosmological spacetimes aimed at both cosmologists and numerical relativists. Second we survey the present body of work focusing on general relativistic simulations without approximations organised according to the cosmological history  from cosmogenesis through the early hot Big Bang to the latetime evolution of universe. In both cases we discuss the present stateoftheart and suggest directions in which future work can be fruitfully pursued.
['General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
Cosmology using numerical relativity
{'Cosmology using numerical relativity': "Josu C. Aurrekoetxea 1 , Katy Clough 2 , Eugene A. Lim 3  \n1 Astrophysics, Denys Wilkinson Building, University of Oxford, Keble Road, Oxford OX1 3RH, United Kingdom.  \n2 School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E1 4NS, United Kingdom.  \n3 Theoretical Particle Physics and Cosmology Group, Physics Department, King's College London, Strand, London WC2R 2LS, United Kingdom.  \nContributing authors: josu.aurrekoetxea@physics.ox.ac.uk; k.clough@qmul.ac.uk; eugene.lim@kcl.ac.uk;", 'Abstract': 'This review is an up-to-date account of the use of numerical relativity to study dynamical, stronggravity environments in a cosmological context. First, we provide a gentle introduction into the use of numerical relativity in solving cosmological spacetimes, aimed at both cosmologists and numerical relativists. Second, we survey the present body of work, focusing on general relativistic simulations without approximations, organised according to the cosmological history - from cosmogenesis, through the early hot Big Bang, to the late-time evolution of universe. In both cases, we discuss the present state-of-the-art, and suggest directions in which future work can be fruitfully pursued.', '1 Introduction': "Einstein's theory of General Relativity (GR) revolutionised our understanding of cosmological evolution, and has provided the basis for modern precision cosmology. However, most studies rely heavily on the symmetry provided by large-scale spatial isotropy and homogeneity, as encapsulated by the spatially flat Robertson-Walker (RW) metric,  \nds 2 = -dt 2 + a 2 ( t )( dx 2 + dy 2 + dz 2 ) . (1)  \nWhere such a description holds, a preferred frame exists, and we can talk meaningfully about cosmological time and space , as defined by some observers aligned with the symmetry. The Hubble parameter H ( t ) = ˙ a/a , describing the curvature of the metric in the time direction, sets the length scale of the problem at some given time t (as L ∼ H -1 ). The symmetry in the metric relies on assumptions about both the underlying curvature and topology of the spacetime, and the distribution of matter. In particular, present day observations have shown that matter in the universe must be roughly spatially homogeneously distributed, with no O (1) overdensities on Hubble scales. Since, at the largest scales at least, this description does appear to correspond well to the current observable universe (see e.g. [1]), most cosmological simulations can be accurately approximated by splitting the background evolution of the scale factor a ( t ) from the dynamics of the matter fields. For horizon-scale modes where the density is low and the fields are linear, a perturbative approach can be used, which is the usual approach for inflationary spacetimes (see e.g. [2]). For low density but highly non-linear field dynamics at the Hubble scale, one can numerically evolve the matter equations of motion on an expanding RW background, using the spatial average of the matter density to drive the evolution of the scale factor for consistency. If the perturbations are large, but sub-horizon and non-relativistic (i.e. local backreaction matters), then the dynamics can be approximated using Newtonian gravity on an expanding RW background, a strategy that is the leading approach in large-scale N-body simulations of late-time cosmological evolution (see e.g. [3, 4]).  \nWhile numerical relativity (NR) is a standard tool in the computation of accurate gravitational waveforms in compact object astrophysics (i.e. black holes (BHs) and neutron stars (NSs)), its application to cosmology is not yet widespread. Recent works that implement improvements in gauge choices, initial conditions and stability have nevertheless been pushing forward our understanding of the universe, building on a number of pioneering early explorations in the field. Cosmologists are now beginning to harness the considerable power of NR to explore strong gravity problems that are beyond the reach of standard perturbative treatments.  \nGiven the success of models relying on background RW symmetries, and the additional cost and complexity involved in NR, one should question the need to go beyond them to fully general relativistic simulations with no preferred frame. In general, three conditions should be met for an NR simulation to be justified:  \n- 1. Local backreaction is large with relativistic matter exhibiting strong inhomogeneities in the density δρ/ρ average ∼ 1.\n- 3. Spacetime is dynamical and absent of symmetries that could simplify the description.\n- 2. The inhomogeneities are coherent on Hubble scales L ∼ H -1 or over Hubble times T ∼ H -1 .  \nThe relation of cosmological NR simulations to other regimes where approximations can be made is illustrated in Fig 1.  \nOne important but less obvious example is work that seeks to rule out the accumulation of secular or non-linear effects that may be lost when adopting a perturbative description. An example of this is work to address concerns about so-called backreaction effects being unaccounted for in cosmological evolutions. As we will describe later in this review, these seem largely unfounded, with effects at percentage levels when studied systematically in full GR. Nevertheless, the simulations helped to advance discussions about how to correctly compare the non-linear and perturbative results, and the impact on next generation surveys and their systematics. Other examples of problems requiring NR that we will see in this review include the study of cosmological singularities, mechanisms like inflation and slow contraction that aim to act as smoothers of initially inhomogeneous spacetimes, and non-perturbative dynamics in reheating or primordial black hole (PBH) formation.  \nIn this review, we focus on fully general relativistic simulations without approximations (i.e. not assuming a certain background or perturbative description). We provide useful background in Sec. 2. Recognising that most cosmologists are not familiar with the methods of numerical relativists, and vice versa , we start by providing two reviews of the topic from different angles - first we explain cosmology  \nFig. 1 : In cases in which overdensities become large δρ/ρ average ∼ 1 on the length scale set by Hubble L ∼ H -1 , and in the absence of any obvious symmetries, NR methods are required to make progress. BG and AF stand for background and asymptotically flat , respectively.  \n<!-- image -->  \nto numerical relativists, and then we explain NR to cosmologists. We then review common technical and conceptual challenges of using NR in cosmological applications, which are themselves open topics of research. The resolution of some of these challenges may lead to new techniques that can be used in other applications of NR, and also to deeper insights about the nature of the cosmological spacetime in which we exist.  \nNext we provide a review of past work, divided into three broad sections in order of time, as illustrated in Fig. 2. The first, Sec. 3, covers the most distant past of the universe, and questions about how it may have begun, including work on cosmological singularities, and attempts to reconcile our homogeneous universe with an inhomogeneous beginning. The second, Sec. 4, describes strong gravity phenomena that may have occurred when our universe was already approximately homogeneous and isotropic on large scales, such as the non-perturbative dynamics in preheating and PBH formation. Finally Sec. 5 covers the late universe, and in particular quantification of the potential effects of full GR on results from the approximate descriptions that are used in most large-scale cosmological simulations. We close with a summary table of the work done to date and a list of open directions.  \nWe would highlight several excellent reviews on related topics, which complement the overview provided here, including [5-13].", '2 Background on numerical relativity and cosmology': 'We start by providing sufficient background of both NR and cosmology to frame the later discussions on specific works. This is by no means a comprehensive description of either NR or cosmology, but aims to highlight the key features of each field that are relevant in the context of the other. With this background in place, we will discuss the key technical and conceptual issues associated with cosmological spacetimes in strong gravity regimes, namely boundary conditions, initial data, gauge choice, interpretation and units.', '2.1 Cosmology for numerical relativists': 'We provide here only a brief description of cosmology in the context of NR. A good starting point for those new to cosmology are textbooks by Baumann [2], Weinberg [14] and Mukhanov [15].', 'Revelation: The late universe': "- 1. What is the impact of including GR effects on systematics for next generation cosmological surveys?\n- 2. Does the assumption of periodicity that is used in most simulations impact on the (lack of) development of non-FLRW features?\n- 3. Are there other GR effects that could explain cosmological tensions?\n- 4. What are the signatures of modified gravity effects in the large scale structure statistics, e.g. non-linear screening?  \nAcknowledgements. We acknowledge helpful input from many of the authors of the papers presented here. JCA acknowledges funding from the Beecroft Trust and The Queen's College via an extraordinary Junior Research Fellowship (eJRF). KC acknowledges funding from from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 693024), an STFC Ernest Rutherford Fellowship project reference ST/V003240/1 and STFC Research Grant ST/X000931/1 (Astronomy at Queen Mary 2023-2026). EAL acknowledges support from a Leverhulme Trust Research Project Grant.", '2.2 Numerical relativity for cosmologists': 'As for cosmology, we provide here only a brief description of NR in the context of the work that we will review. A good starting point for those new to NR is the student-friendly text by Baumgarte & Shapiro [28] and for those looking to understand the concepts here in more depth the textbooks by Alcubierre [29], Baumgarte & Shapiro [30], Gourgoulhon [31], and Shibata [32].  \nThe basic principle of NR is to solve the coupled matter-metric system of GR as a Cauchy or initial value problem, whereby we specify initial data for the matter and the metric on some space-like hypersurface, and then evolve it forward in time (defined according to some chosen observers); that is, we specify the initial scenario, and then use the computer to find out what happens next. Such a formulation is naturally provided by the ADM decomposition of the Einstein equations, which was developed in the 1960s by Arnowitt, Deser, and Misner [27].  \nA standard question when confronted with any formulation of PDEs is whether they are well-posed , which is the demand that a unique solution exists that depends continuously on the initial data. In the perturbative cosmology of the previous Sec. 2.1, one separates the background evolution and the perturbation evolution equations, both of which are well-posed. In fully non-linear NR, the well-posedness of the formulation is less obvious and problems with well-posedness in NR have a long history. In particular, one of the early obstacles in NR was that the ADM formulation [27] was ill-posed, rendering its direct numerical implementation unstable. Fortunately, this problem can be surmounted by reformulating the evolution equations into one of several possible (locally) well-posed forms. For a detailed introduction to well-posedness in NR, see the review article by Hilditch [33]. One popular formulation is the so-called Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism [34, 35] based also on an earlier work with Nakamura, Oohara and Kojima [36]. An alternative formulation is the Generalised Harmonic Coordinates (GHC) formalism [37], which is well-known for its use in successfully evolving binary BH mergers for the first time by Pretorius [38], but was applied to studies of the cosmological singularity before this by Garfinkle [39]. Another alternative is the tetrad formulation developed by Uggla et al. [40], used by Garfinkle for the study of cosmological singularities [41] and recent works on bouncing cosmologies (see Sec. 3.3), which allows scaling out of the cosmological scale H . This helps the numerics to follow longer timescales during which H is changing over a wide range of scales, as at the end of inflation or contraction  \ntime', 'Gauge in NR': '<!-- image -->  \n<!-- image -->  \ndensity space space Fig. 3 : Gauge in cosmology versus NR. (Left:) For cosmologists a change in gauge implies a small perturbation to the slicing, but with the background FLRW slice about which the perturbation is defined always aligned with the symmetry of the spacetime. This can transform a constant density slice into one in which small perturbations are observed. (Right:) In the NR case, a change of gauge is a fundamentally different slicing of the spacetime, which therefore potentially results in a (non-perturbative) change in the density observed on each constant-time slice. It is important to note that when there is no underlying symmetry, the notion of a preferred slicing is lost, and thus there is no choice of gauge in which we are able to recover FLRW-like observers.  \n<!-- image -->  \ntowards a singularity, but necessitates the solution of an elliptic equation for the lapse at each step. There are some other early works that formulate hyperbolic equations with tetrads [42, 43], which were mostly abandoned when coordinate basis formulations achieved success with binary BH mergers, but may be worth revisiting to address some of the challenges of cosmology.  \nIn this section we will mainly discuss NR in the context of the BSSN formalism, but the principles are similar for most formulations. The idea of slicing up spacetime into space-like hypersurfaces and then evolving them is familiar to most cosmologists, where the RW metric imposes a fixed foliation and coordinates upon each hyperslice as discussed in Sec. 2.1. In NR, the choice of foliation and the coordinates imposed on it are usually dynamical - they respond to the local evolution of the matter and metric. This is a necessity for a well-posed evolution, but can also help with avoiding singularities that form and tracking changing scales over time on a finite numerical grid.  \nWe now briefly review the ADM formalism [27] and the way in which it is adapted for NR. The ADM formalism projects the 4D metric g µν into a foliation of spatial hypersurfaces with metric γ ij as in Fig. 4, with surfaces labelled by a constant time coordinate t . The normal vector to the slices is n µ = (1 /α, β i /α ). The coordinates x i runs over the three spatial directions, and the metric expressed in a coordinate basis adapted to the slicing is  \ng µν dx µ dx ν = -( α 2 -β i β i ) dt 2 +2 β i dx i dt + γ ij dx i dx j . (20)  \nIn NR, the coordinates t and x i refer to the numerical time and grid locations respectively, and, whilst they can be related to particular physical observers, these are not in any sense preferred in the same way as cosmological observers, or the asymptotic observers in a Schwarzschild spacetime. The metric degrees of freedom describing proper lengths on the hypersurface are encoded in the 3-metric γ ij , while α and β i are the lapse and shift variables respectively, which are generally thought of as gauge functions that define the choice of foliation and the coordinate labels imposed on it, although the clear separation of  \ndensity  \nFig. 4 : Foliation of a 4D spacetime into a 3D spatial hyperslice, and a time-like normal to that slice. The gauge variables - the lapse α and shift β i - are also illustrated. In this picture space is represented as a two dimensional surface, whereas in full GR each spatial slice is a 3D volume.  \n<!-- image -->  \nthe metric degrees of freedom into physical and gauge is not easily defined. Note that we have implicitly chosen the coordinate basis defined by the Cartesian grid x µ which is common in NR but not the only choice (for example the tetrad basis mentioned above or the dual foliation approach of Hilditch [44]).  \nGiven that the equations of motion for the metric are second order in time, in NR we introduce the time derivative of the metric, in the form of the extrinsic curvature K ij , to reorganise the equations as a coupled, first order in time system. The extrinsic curvature is equivalently defined as the Lie derivative of the metric along the normal to the slice,  \nK µν = -1 2 L n g µν , (21)  \nor as the projection of the quantity ∇ µ n ν into the spatial slice K µν = γ σ µ γ ρ ν ∇ σ n ρ . The spatial components of the extrinsic curvature in the adapted basis  \nK ij = 1 2 α ( D i β j + D j β i -∂ t γ ij ) , (22)  \ntogether with γ ij , forms the phase space of metric degrees of freedom.  \nPhysically, K ij defines how curved the 3D spatial slice is within the 4D spacetime in which it is embedded. A simple example of the difference between intrinsic curvature and extrinsic curvature is the surface of a cylinder - the 2D surface is flat; were we to draw triangles on it, the angles would add up to 180 degrees. However, we consider it to be curved because the normal to the 2D surface changes direction in the 3D space in which it is embedded - this latter concept of curvature, related to how the direction of the normal vector changes direction along the surface (recall that K µν ∼ ∇ µ n ν ), is the extrinsic one, see Fig. 5. Again cosmology provides a nice simple illustration of the concept, as for the RW metric Eqn. (1) we can identify the extrinsic curvature as being related to the time derivative of the spatial metric, and so to a ( t ) and H ( t ), as  \nγ ij = a 2 δ ij , K ij ∼ -Hγ ij , (23)  \nand the trace K = γ ij K ij = -3 H . Note that some groups use the opposite sign convention (here we use that positive and negative K corresponds to collapse and expansion, respectively). This reminds us that although our universe is roughly spatially flat, it is not at all flat in 4D spacetime, but actually curved in the time direction, by an amount quantified by H . We perceive the curvature of the metric in the time direction as an expansion of space in time, but it is really that geodesic observers are moving on a strongly curved spacetime and thus experience geodesic deviation along their trajectories.  \nThe definitions so far are only geometric, and do not encode the physical dynamics of the Einstein equations. To do so the equations themselves must be decomposed into space-like and time-like parts.  \nFig. 5 : The extrinsic curvature K µν encodes the time derivative of the spatial metric γ ij , which can be related to how the direction of the normal vector changes along the surface, projected into the surface K µν = γ σ µ γ ρ ν ∇ σ n ρ . (Note that n µ is orthogonal to the hypersurface in the Lorentzian sense since this is a spacetime diagram).  \n<!-- image -->  \nDoing this we obtain an evolution equation for the extrinsic curvature K ij  \n∂ t K ij = β k ∂ k K ij + K ki ∂ j β k + K kj ∂ i β k -D i D j α + α ( (3) R ij + KK ij -2 K ik K k j ) +4 πα ( γ ij ( S -ρ ) -2 S ij ) (24)  \nas well as the Hamiltonian and three momentum constraints  \nH ≡ R + K 2 -K ij K ij -16 πρ = 0 , (25)  \nM i ≡ D j ( γ ij K -K ij ) -8 πS i = 0 , (26)  \nwhich must be satisfied on every hypersurface. The decomposed components of the matter stress-energy tensor (those measured by normal observers) are defined as  \nρ = n µ n ν T µν , S i = -γ iµ n ν T µν , S ij = γ iµ γ jν T µν , S = γ ij S ij . (27)  \nThe equations of motion for the matter fields (and for a fluid, their equation of state), close the system of equations.  \nAgain comparison of these ADM equations to the FLRW case is instructive. We see that the momentum constraints are trivially satisfied via isotropy, and the Hamiltonian constraint is simply the first of the Friedmann equations in Eqn. (3). Whilst this is often thought of as an evolution equation for the scale factor, it is in fact a constraint, and the true dynamical equation for the metric quantity a is the second Friedmann equation of Eqn. (3) as discussed in Sec. 2.1 4 . The more complicated evolution equation for K ij , Eqn. (24), is an extension of the 2nd Friedmann equation to the case where the quantities are spatially inhomogeneous and the gauge is non-trivial. For example, if K xx is greater than K yy this implies that (according to our chosen observers) the spacetime in the x direction is collapsing faster than in the y direction.  \nWith the degrees of freedom and their equations of motion defined, the system can in principle be discretised and solved on a computer using standard numerical integration and differentiation techniques (see e.g. the relevant chapters of [29, 30]). However, as discussed above, due to the additional complication of ensuring the equations are well-posed, the ADM equations are modified and auxilliary variables may be defined depending on the formulation, see Hilditch [33] for a review.  \nThe basic conceptual process of setting up and running a NR simulation is as follows, and we will provide further discussion of each step in the following sections:  \n- 1. Define the spatial boundary conditions, which may include assuming a certain topology of our universe e.g. T 3 × R - see Sec. 2.3.  \n- 2. Specify initial conditions that both satisfy the Hamiltonian and momentum constraints of GR, and are physically the correct ones for the system to be studied - see Sec. 2.4.\n- 3. Define a gauge/coordinate evolution that is stable and well adapted to the problem - see Sec. 2.5.\n- 4. Extract and interpret meaningful diagnostic quantities and observables - see Sec. 2.6.  \nA key difficulty in any numerical simulation is the resolution of different length scales, since computers have a finite precision with which they can represent numbers. For example, if we want to study small length scale effects over long timescales (relative to the light crossing time of the physical length scale), it will be challenging. This is both a feature and a bug, as we have already seen that NR is most needed where the length scales of interest are similar and therefore perturbative approaches fail. However, there will still be many cases where the length scales in the problem change significantly in space or time, and therefore the numerics can break down. In particular, cosmological spacetimes with local overdensities introduce an additional global scale (i.e. the Hubble parameter H ) into the numerical system. This creates a challenge when our goal is to study local physics with a different scale. In many cases we expect strong backreaction to be important on small scales, and while it is possible that the dynamics can decouple from the background cosmological evolution, there are intermediate cases where both scales need to be fully resolved. In cases where we need to follow a range of scales in different regions, one can use Adaptive Mesh Refinement (AMR) techniques to gain one or two orders of magnitude, but it is even better to use clever choices of the evolution variables that scale out expansion, as in the tetrad formulation of Uggla et al. [40] and BSSN-like formulations as used by Zilhao et al. [45]. See also Daverio et al. [46] for a proposed scheme for cosmological spacetimes that eliminates the fluid components through the constraint equations.', '2.3 Boundary conditions': 'The most important difference between NR applications to astrophysical systems and to that of cosmological systems is that, at the largest scales, the universe is not asymptotically Minkowski/flat. For most of its evolution, it is well described by an FLRW universe with a characteristic scale at time t set by the Hubble parameter H ( t ). The connection between the astrophysical use of NR and the reality of our cosmology is that the length and timescales for most astrophysical phenomena are much shorter than the Hubble length and time ∼ H -1 , and so they do not feel the cosmological evolution. One can then safely assume that spacetime is asymptotically flat, for example, in a local phenomenon such as a BBH merger. On the other hand, in an NR simulation of cosmological spacetimes, we are often interested in the global structure over cosmological timescales, e.g. how are the matter fields distributed? Is the spacetime inflating? etc. Therefore, the fact that the universe is not Minkowski on average will play a role and the choice of spatial boundary may impact on the physical result of the simulation.  \nThe most common and straight-forward solution is to apply periodic boundary conditions, usually with the idea of treating the computational domain as being a typical patch of the entire larger universe, rather than being truly topologically a torus on the relevant scales, see Fig. 6. In late-time post-Big Bang Nucleosynthesis cosmological evolution this is a well-justified choice as the spacetime is observed to be statistically homogeneous and isotropic 5 . In the early universe the choice is less easy to justify we simply do not know what the topology of our universe is on scales beyond the cosmological horizon. Imposing periodic boundary conditions in Cartesian coordinates restricts the topology of space to be a 3-Torus, and hence the Ricci scalar of the spatial hypersurface (3) R has to be negative somewhere or (3) R = 0 everywhere if there are no BHs [47]. Operationally, periodic boundary conditions also impose restrictions on the choice of initial conditions [48] - we will discuss this in Sec. 2.4 later. Finally, periodic boundary conditions impose a discrete symmetry in the size of the domain with a physical scale that is potentially evolving since it depends on the gauge. Thus one must ensure that the scale of the physics being investigated is smaller than the scale of this discrete symmetry, or accept that the results may depend on the choice of boundaries. This issue has recently been raised in the context of studies of inflation in [49], and more generally in N-body studies of cosmological evolution in [50].  \nAlternative approaches have been taken, in which one instead matches to a particular FLRW spacetime asymptotically, see for example [51-53]. There is unfortunately no guarantee that such a choice will be consistent - imagine we choose to match to an asymptotically de Sitter region during a simulation of inflation [54], and some dynamics of the field at the interior of the domain propagate waves out to the boundary. This should disturb the boundary in some way and so rigid de Sitter becomes inconsistent and creates unphysical reflections of the waves. One may be able to fix such a case by applying de Sitter  \nFig. 6 : Periodic boundary conditions are often used in cosmological simulations, with the idea of treating the computational domain as being a typical patch of the entire larger universe, rather than being truly topologically a torus on the relevant scales.  \n<!-- image -->  \nplus outgoing waves as is done for Sommerfeld boundary conditions in an asymptotically flat spacetime. However, the problem is more general - the choice of boundaries has a physical impact if the timescales of interest allow them to communicate with the physical region of interest, and therefore their choice is part of the specification of the physical situation under study.  \nThis is an area in which further work is needed, and new ideas would be welcome. Some options might be to revisit alternative boundary choices that have previously been explored in NR, such as absorbing or constraint preserving boundary conditions [55-60], to try to compactify the domain as in recent work for gravitational wave extraction [61-64], or to impose that the boundary tracks a cosmological horizon or restricts the domain to that bounded by the null lightcones emanating from the boundaries.', '2.4 Initial data for cosmological simulations': 'One of the key difficulties in NR is how to set initial conditions that both satisfy the constraint equations and represent the physical spacetime of interest. In this section, we will focus on the issues that arise in cosmological NR simulations. We refer the reader to the standard NR texts [28-30, 32] for more details on the methodologies for finding initial data in NR.  \nInitial data for the gravitational degrees of freedom consists of the 6 components of the spatial metric γ ij and the 6 components of the extrinsic curvature K ij at each point on the initial hypersurface. For a given initial matter configuration, { γ ij , K ij } must satisfy the Hamiltonian and momentum constraints Eqns. (25-26), which, in the most general case, represent a set of four coupled, elliptic PDEs. Configurations that satisfy these equations are said to be on the constraint surface . The four constraints remove 4 degrees of freedom from the initial data - the remaining 8 (4 of which are physical, and 4 gauge) together with the 4 gauge variables { α, β i } must be chosen according to physical principles or knowledge about the system. A key challenge is that, even when we have an initial solution which is on the constraint surface, it might not reflect the physical problem that we are trying to investigate . This is because, for a given matter configuration, described by the stress-tensor T µν , there may be multiple (and potentially an infinite number of) metric configurations that satisfy the Einstein equations that are not physically equivalent. For example, one can usually add some additional transverse traceless component to the extrinsic curvature K ij (roughly additional gravitational wave content). This will change the physical spacetime, but the constraints can still be satisfied.  \nThere is no one-size-fits-all method to find initial conditions in NR. Nevertheless, there are two standard methods for finding initial conditions - the conformal traverse-traceless (CTT) and the conformal thin sandwich (CTS) approaches [65, 66], which as their names imply, perform a conformal transformation  \nof the 3-metric  \nγ ij = ψ 4 ¯ γ ij . (28)  \nThe simplifying assumption of conformal flatness ¯ γ ij = δ ij immediately removes 5 degrees of freedom from the 6 of the metric γ ij , leaving the the conformal factor ψ to be specified 6 . The 6 degrees of freedom in K ij are decomposed into  \nK ij = A ij + 1 3 γ ij K, (29)  \nwhere A ij can be further decomposed into a traceless-traverse component ¯ A ij TT and a longitudinal component ¯ A ij L ,  \n¯ A ij = ψ 2 A ij , ¯ A ij = ¯ A ij TT + ¯ A ij L . (30)  \nThe 6 components - 1 in K , 2 in ¯ A ij TT and 3 in ¯ A ij L - combined with ψ , represent the remaining 7 components that can be set. Given the 4 constraints in Eqn. (26), this means that we still need to choose values for 3 components. The choice of which components to specify, and which are resolved using the constraints is dependent on both the physical problem at hand, and the simplicity (or not) of finding solutions.  \nIn the special case when the stress-tensor has zero momentum density S i = 0, the momentum constraint Eqn. (26) becomes  \nM i = D j ( γ ij K -K ij ) = 0 , (31)  \nand hence D i K = D j K ij . A simple and often used ansatz is to set K = constant and ¯ A ij = 0 so that the momentum constraint is trivially satisfied. Such a choice leaves ψ to be determined, and this is obtained by solving the Hamiltonian constraint, an elliptic equation in ψ that can be solved with standard multigrid techniques,  \n8 ¯ D 2 ψ -ψ ¯ R -2 3 ψ 5 K 2 + ψ -7 ¯ A ij ¯ A ij = -16 πψ 5 ρ , (32)  \n(where we have restored ¯ A ij and a general conformal metric for completeness). Here ¯ D is the covariant derivative associated with the conformal 3-metric ¯ γ ij and ¯ R is its Ricci scalar.  \nIn more general cases, one in principle has to solve the coupled set of constraints Eqn. (26). Both the CTT and the CTS methods that were initially developed for astrophysical simulations such as BH and NS mergers have been successfully used in cosmological NR simulations. In the former, the freely specified variables are { ¯ γ ij , K, ¯ A ij TT } together with the gauge variables α and β i , while in the latter the freely specified variables are { ¯ γ ij , ∂ t γ ij , K, ∂ t K } and α and β i are the determined variables.  \nIn both of these methods, the conformal factor ψ is determined by solving the Hamiltonian constraint, an elliptic equation Eqn. (32). However, solutions to elliptic equations do not always exist, and when they do, they are not necessarily unique. Indeed, if periodic boundary conditions are imposed as is common in cosmological simulations (see Sec. 2.3), elliptic equations of the form of Eqn. (32) do not possess unique solutions or are necessarily trivial because the Laplacian operator ¯ D 2 is not invertible on a conformally flat spatial slice with periodic boundary conditions. For periodic boundary conditions, a set of integrability conditions must be obeyed such that not only are the constraints in Eqn. (26) satisfied, but their integrals over the periodic domain are zero (more specifically, one needs to ensure that the source of the elliptic equation is in the adjoint of the Laplacian operator, see [67, 68]). For example, for the case where ¯ A ij = 0 and ¯ γ ij = δ ij , the integrability conditions  \n∫ dV ( 2 3 ψ 5 K 2 -16 ψ 5 ρ ) = 0 , ∫ dV ( 2 3 ψ 6 ∂ i K +8 πψ 6 S i ) = 0 , (33)  \nmust be satisfied. This can lead to difficulties in constructing initial conditions for cosmological spacetimes, and limit the choice of initial conditions that can be explored.  \nThe presence of fundamental matter fields plays a crucial role in many cosmological simulations. In such cases, matter is represented by one or more fundamental fields such as the scalar inflaton, which presents a special challenge in setting initial conditions. This is because unlike fluids, which can be represented by macroscopic variables such as energy density ρ and pressure p , these variables themselves  \nare functions of the fields and the metric. For example, for a single scalar field ϕ with potential V ( ϕ ), the energy density and momentum are  \nρ = 1 2 ψ -4 γ ij ∂ i ϕ∂ j ϕ + 1 2 Π 2 + V ( ϕ ) , S i = -Π ∂ i ϕ , (34)  \nwhere the conjugate momentum Π ≡ α -1 ( ∂ t ϕ + β i ∂ i ϕ ). This is a complication because standard methods to ensure the elliptic equation Eqn. (32) has unique solutions, that involve scaling the energy and momentum densities [69, 70], do not easily extend to fundamental fields, especially in the case of a non-trivial potential.  \nThe problem arises due to the maximal principle [71]. Setting ¯ A ij = 0 and ¯ γ ij = δ ij for simplicity and linearizing Eqn. (32) via ψ → ψ 0 + u , we get the following inhomogeneous linear PDE  \n( ∂ j ∂ j -q ( x ) ) u = σ ( ψ 0 ) , q ( x ) = 5 12 ψ 4 0 K 2 -10 πψ 4 0 ρ , (35)  \nwhere the source σ ( ψ 0 ) = -2 ψ 5 0 ρ -∂ j ∂ j ψ 0 + (1 / 12) ψ 5 0 K 2 . If the boundary is Dirichlet and q ( x ) > 0 everywhere in the domain and its boundary, then the maximal principle implies that the only solution is u = 0 and hence a unique solution for ψ in Eqn. (32) exists. Conversely, if q ( x ) < 0 somewhere, then this cannot be guaranteed and a unique solution may not exist. A standard trick to render the solution unique is to rescale ρ → ˜ ρψ -5 , which switches the sign of the ρ term in σ and ensures that q ( x ) > 0. If ρ itself depends on ψ and ϕ , such a rescaling may not always allow ϕ to be simply reconstructed, although in Corman & East [72] a rescaling of the field itself was successfully employed for a massive potential, along with bespoke analytic solutions that ensure the momentum constraint is always satisfied as proposed by Garfinkle & Mead [67].  \nIf rescaling is not possible, one can still attempt to numerically solve the constraints to find a solution (if it exists) - if a trial solution is close to a solution and the perturbations are weak, this method can still work. Alternatively, in Aurrekoetxea et. al. [48], it was proposed that instead of specifying K and solving the Hamiltonian constraint Eqn. (32) for ψ as an elliptic equation, one instead specifies ψ and then solves for K as an algebraic equation (up to a sign). This so-called CTTK method has been successfully employed in cosmological NR problems [73-76] and was inspired by the Friedmann equations where H ∼ K varies with ρ . It also works for BBH scenarios [77, 78], but note that as a result of the spatial variation in K , solving the momentum constraint cannot be avoided.  \nTo date most works use a conformally flat metric ¯ γ ij = δ ij as an initial condition. Whilst in the context of simulating cosmological spacetimes this may seem a natural choice, we should be mindful of several things. As pointed out by Garfinkle & Mead [67], this imposes a constraint on the choice of A ij and S i - there exist matter configurations that do not source a conformally flat spatial metric hyperslice (in particular, those with significant anisotropy). One good example is a uniform momentum density in a single direction, which does not obviously conflict with the periodicity imposed but is prevented by the integrability conditions in the conformally flat case. This is a strong simplification in terms of the possible initial conditions and may be important in simulations to test the robustness of inflation and bouncing spacetimes (see Sec. 3.2 and Sec. 3.3). Going beyond conformal flatness requires a generalisation of the methods described above, but should not be conceptually problematic - this would be an interesting avenue for future work and permit studies of more general cases in cosmological spacetimes.  \nOther tricks have been employed to solve the Hamiltonian and momentum constraints by making judicious choices of the matter fields, such as solving for the momentum of the scalar field or its the energy density given an inhomogeneous choice for the conformal factor ψ (e.g. [79, 80]), or introducing extra fields to cancel out or simplify the overall matter sources in the constraints (e.g. [81-83]). Such methods are a good choice where we are not too concerned about the exact form of any perturbations we are introducing, and just want constraint satisfying data in which perturbations are present in the metric.  \nMore work on connecting the physical properties of the initial spacetime to the choices of free variables would be valuable for cases where we would like more control over the data. We also note that at present there seems to be no fully general, publicly available initial condition solver for cosmological spacetimes, unlike codes for asymptotically flat spacetimes of which there are several 7 . One exception is FLRWSolver , developed by Macpherson et al. , but this is specialised to initialise data for cosmological perturbations arising from inflation for studies of late-time cosmology. More open-sourcing of such tools would lower the barrier for new researchers undertaking NR simulations, which would ultimately benefit the field.', '2.5 Gauge choices': 'There are two key considerations that set the requirements for a choice of gauge in numerical simulations - firstly, and most importantly, is the stability of the evolution. The second consideration is that we would like to efficiently sample the spacetime, and want to distribute the coordinates across the region of interest in a way that captures the relevant physical scales, but does not oversample them. As discussed above, in a BSSN-type formalism, the gauge freedom is broadly encapsulated in the lapse and shift functions, α and β i .  \nIn the case of zero shift β i , the lapse function α determines the ratio between the proper time experienced by a normal/Eulerian observer and the coordinate time that passes for the observers at fixed spatial coordinates. A comparison of the ADM metric Eqn. (20) and RW metric Eqn. (1) reminds us that in the spatially homogeneous case, the lapse is related to the reparameterisation of the time coordinate, for example, moving from conformal time η to proper time t as dt = adη , or simply choosing the unit of time. In an inhomogeneous spacetime, the effect is more complex and in particular a spatially varying lapse results in an acceleration of the coordinate observers relative to the normal ones.  \nThe shift vector β i determines how the coordinate points relating to each observer on the slice are relabelled as we move forward in time as illustrated in Fig. 4. Again in a homogeneous case this is easy to visualise: imagine β x = -1 and so we relabel our x coordinate as x + dt after a time dt elapses. This is just a Galilean transformation in the x coordinate, and would make an object that was originally at a fixed coordinate appear to translate with a coordinate speed of dx/dt = 1. However, the spatially varying shift case is less easy to interpret and introduces a physical distortion of the spacetime - we see from Eqn. (22) that its gradients contribute to the extrinsic curvature K ij , for example. In most NR simulations both α and β i are functions of both x i and t .  \nSetting α = 1 and β i = 0 is a possible choice of gauge - this gives a geodesic slicing in which the numerical grid corresponds directly to the normal or geodesic observers. In the context of NR simulations of cosmology, if deviations from an FLRW spacetime are small, then this gauge is equivalent to setting the synchronous gauge in perturbation theory (see Sec. 2.1). Unfortunately, despite its simplicity, geodesic slicing is a poor choice of gauge when strong localised overdensities and/or BHs form, since gravitating systems tend to focus geodesic congruences, causing coordinate singularities. Therefore, both the gauge is usually allowed to vary spatially and in time during a simulation. In particular, in the moving puncture approach, evolution equations are imposed on ∂ t α and ∂ t β i as functions of the metric variables, or equivalently in the approach of GHC, some appropriate gauge source functions H µ are specified. More details on these approaches are given below.  \nThe choice of gauge evolution is usually tied to the physical problem at hand, and we will discuss particular NR gauge choices used in cosmological works in the remainder of this section. In Sec. 2.6 we will discuss the difficulties in interpretation that the use of a dynamical gauge introduces. However, let us emphasise already here that where the spacetime is highly inhomogeneous, the notion of a preferred time-like observer has already been lost. Thus even in the case where we chose the lapse and shift to be trivial, there is no right gauge and we would not be able to recover FLRW observers, as we discussed in Sec. 2.1, see also Fig. 3.  \nThere is no magic bullet in the choice of gauge in NR, but there are now some well-established principles for what works. A good choice of gauge should typically aim to do some or all of the following (a) avoid hitting singularities in a finite time (b) maintain well-posedness of the system of PDEs (c) settle into an approximately stationary state when the spacetime possesses an approximate time-like Killing vector (d) avoid stretching the hyperslice by driving the spatial metric towards conformal flatness (e) ensure the discrete timesteps respect the CFL condition for numerical stability - i.e. ensure that the ratio of proper time and proper distance is less than the O (1) number determined by the time integration method used dτ/dl < O (1). Choosing a gauge in NR is often a matter of trial and error (and some good luck).  \nCosmological spacetimes that are uniformly expanding are often more forgiving than BH ones when it comes to the choice of gauge - as long as overdensities remain small, geodesic slicing, with the lapse α = 1 and the shift β i = 0 can be stable, and this choice has been employed successfully in several works (e.g. [79, 84]). As discussed above, in the limit of a homogeneous spacetime, this choice matches the standard FLRW synchronous gauge. The stability arises because the expansion of the normal observers counteracts the tendency for them to focus as a result of overdensities. However, once significant overdensities form, focusing of geodesics will lead to coordinate singularities forming in finite time and one of the standard NR gauges or extensions thereof must be used.  \nThe standard approach to evade these singularities is the moving puncture approach, where the lapse and shift variables are given dynamical evolution equations. For the lapse evolution, most simulations use a particular numerically robust class of foliation defined by the so-called Bona-Masso slicing condition [85]  \n∂ t α = -α 2 f ( α ) K , (36)  \nwhere K = Tr( K ij ) is the trace of the extrinsic curvature and f ( α ) is some (normally positive definite) function of α . For K > 0 (contracting regions that normally occur around overdensities), α → 0 is driven to small values to prevent the coordinate observes from focussing and creating a coordinate singularity (or reaching a physical one in the case of BHs). Often this collapse of the lapse is thought of as slowing the passage of time so one never reaches the singularity, but its effect is mainly to create a gradient in the lapse that generates an outward acceleration of the normal observers relative to the overdensity, thus stabilising them relative to freely-falling geodesic observers. One can also think of it as an attempt to dynamically drive K to zero.  \nSimilarly, imposing a specific evolution for the shift β i can help to avoid singularities and stabilise the evolution - one effective condition is the Gamma-driver [86, 87]  \n∂ t β i = 3 4 B i , ∂ t B i = ∂ t ˜ Γ i -ηB i , (37)  \nwhere η > 0 is a dimensionful constant of order ∼ 1 /M that damps the shift evolution, and B i (not to be confused with the B i variable in Eqn. (7)) is an auxiliary variable. The Gamma-driver avoids slice stretching by seeking to minimise ˜ Γ i ∼ -∂ j ˜ γ ij , and hence preventing the formation of coordinate singularities. Using Eqn. (36) with f ( α ) = 1 /α (1 + log slicing) and Eqn. (37) combined is known as the moving puncture gauge , which allows for the stable evolution of moving BHs. This gauge permitted the successful simulations of BH mergers [88, 89], without the explicit excision of singularities that is required in the GHC formalism 8 [37, 38]. Both the standard GHC and moving-puncture gauges have been successfully applied to NR evolutions of cosmological spacetimes, including cases in which BHs form [72, 90-93] 9 .  \nOne can encounter problems when applying the BH gauges to cosmological spacetimes. For example, in the moving puncture gauge, because α is driven by sgn( K ) so as to slow down evolution in regions of collapsing spacetime when K > 0, in cosmological expansion K < 0 can lead to the rapid growth of the lapse α ≫ 1. Then, each numerical time-step leads to a large leap in cosmic time, which can cause instability and accuracy issues, especially if the lapse varies significantly across the domain and other timescales of interest exist in the problem. For example, in the context of reheating, the oscillation of the field around the minimum can become under-resolved when the lapse grows too large.  \nTo circumvent this issue, one can implement a more cosmological spacetime-friendly gauge. In Giblin & Tishue [97], a modification of the 1+log slicing condition was used in which  \n∂ t α = -2 α ( K -⟨ K ⟩ ) , (40)  \nwhere ⟨ K ⟩ is the spatial average of K over the hyperslice - the idea is to set the rate of change of α such ∂ t α = 0 when space is homogeneous, since K = ⟨ K ⟩ . This gauge choice approximately ensures that simulation time is identified with the cosmic time coordinate rather than conformal time. A similar gauge was used by Kou et al. [98], given by  \n∂ t α = -2 τα 2 1 + α 2 ( K -⟨ K ⟩ ) (41)  \n□  \nx  \nµ  \n= 0  \n,  \n(38)  \nIn the GHC formalism, the idea is to drive this wave equation with a set of  \nsource functions  \nH  \nµ  \n[37, 39]  \n□ x µ = H µ , (39)  \nwhere H µ is evolved dynamically with a set of differential operators L ( µ ) that can be freely specified as a function of the spacetime coordinates, the metric and its derivatives, and the source functions and their derivatives. Their exact form depends on the physics being investigated and can be chosen in such a way as to achieve slicings similar to those discussed for the moving puncture gauge. Eqn. (39) and its derivatives then form a set of constraint equations, which can be substituted into the Einstein equations to form a set of well-posed equations of motion for the metric.  \nwith β i = 0. This is a version of the radiation gauge g ij Γ ρ ij = 0, and has the late-time asymptotic solution K → ⟨ K ⟩ and ∂ t α = 0 for an approximately homogeneous spacetime. Elley et al. [75] introduced the following gauge obtained by setting f ( α ) = e -α /α in the Bona-Masso condition  \n∂ t α = -αe -α K + β i ∂ i α . (42)  \nAgain, the idea is to prevent a runaway growth of α in homogeneous (and hence expanding) regions, while still avoiding the focusing of observers via the Bona-Masso condition in collapsing regions.  \nAn alternative approach in the tetrad formulation proposed by Uggla et al. [40], requires choosing constant mean curvature slices with homogeneous (but time-dependent) K ( t ) to foliate the spacetime. Such a choice tipically requires solving an elliptic equation for the lapse α at every numerical time step, since the aim is to achieve  \n∂ t K ∼ -D 2 α = 0 , (43)  \nwhich is potentially more computationally expensive and may lead to issues of uniqueness and existence of solutions. Nevertheless, it has the value of possessing a clear physical interpretation - these are observers for which local expansion is homogeneous even if the matter distribution is not (though clearly these will not generically be geodesic observers). Other interesting recent propositions are the use of mean curvature flow, as in [99], or a reference metric approach as in [100]. The former requires implicit timestepping methods, which are less commonly used than explicit ones. However, the additional technical complexity may be worthwhile for achieving better gauge behavior in certain cases.  \nIn general, we view the question of gauges in cosmology as having been insufficiently studied. There is a tendency (including by the authors of this review!) to reuse existing methods for BH spacetimes rather than develop and explore new and potentially more well-adapted paradigms. We hope to see new developments in this area in the future.', '2.6 Interpretation of quantities in an inhomogeneous spacetime': "Given the freedom and dynamical construction of the foliation of spacetime into spatial hypersurfaces described in the previous section, it is clear that any interpretation of the numerical results must take into account the potential dependence on the gauge.  \nOne issue is that for any higher-ranked tensor beyond the spacetime scalar, the component values depend on the basis choice - for example for X µ the physical meaning of X 0 , X 1 etc depends on what basis the index µ is associated with. Even when accounting for that (or using true spacetime scalars [49, 101]), what you see on the numerical grid at some coordinate time t still depends on the choice of foliation. This is particularly important when comparing the numerical results from different simulations. In principle, once a simulation is complete, one can 'refoliate the spacetime' via coordinate transforms into any foliation one might desire, but in practice this is extremely difficult. It may also be the case, given the finite extent of the domain in time and space, that simulations with different foliations capture substantially different regions of the spacetime.  \nAn important distinction is between what we will refer to as gauge dependence and slicing dependence . Consider the specific example of the energy density of the ADM decomposition ρ = n µ n ν T µν where n µ is the unit vector orthonormal to the hypersurfaces of the foliation. This is a spacetime scalar, and therefore it must be gauge independent. However, it is the energy density measured by a particular observer, specifically the one following the normal congruence n µ , and therefore if we change our slicing, and hence the physical trajectory that we label n µ , we change the physical scalar that we calculate. Therefore the value of ρ that we measure depends on the slicing - it is slicing dependent . If we were able to track the original normal observer in the new gauge, and calculate the energy density associated with them, then the value would remain the same, but this is not usually done in practise. Mapping between specific observers in inhomogeneous spacetimes is not straightforward (especially if one changes the initial data), and a simulation in a different gauge will normally be associated with a different set of observers, and their associated observables.  \nA related problem is the tendency to use average quantities across the spatial slices as some measure of the evolution. We must be careful in interpreting these since there is no possible time-like observer who can measure such quantities - only a god-like observer who can freeze time and move around the spatial slice would ever be able to measure a quantity at space-like separated points. Even if the quantity is a spacetime scalar, its integral over the surface is slicing dependent, and so the time evolution may say  \nmore about the gauge than the physics. In particular, seeing some quantity blow up can often be a sign of a coordinate singularity forming, rather than a physical instability. This observation is particularly important in the context of quantifying cosmological observations, which are made on a lightcone and not spatial hypersurfaces.  \nThe most natural and well-motivated physical quantities to track in an evolution are the ones measured by geodesic observers. This still involves some ambiguity - e.g. should they be comoving with the expansion on the initial slice or on the final slice (which would require them to be tracked backwards in time)? Such methods have been employed in late-time cosmological simulations, as we will detail in Sec. 5, where precision is important, but in early universe cases the use of quantities associated with coordinate observers and tracking of spatial averages are the common diagnostics. This is a reasonable approach where the question is fairly binary - e.g. does inflation occur somewhere in the spacetime, or not? Does a BH form (as measured by an apparent horizon forming) or not? However, in cases where the question is more precise, e.g. in which regions that are causally connected to some section of the initial data does inflation occur? , then the diagnostics need to be carefully formulated. The development and use of geodesic-based diagnostic tools in order to facillitate comparison in different gauges is an area where further effort could be usefully expended.  \nA common observable for which accurate predictions are crucial is the gravitational-wave spectrum. At the linear level, a non-local gauge-invariant GW density ρ GW can be constructed by averaging over a region larger than the wavelengths. In the FLRW framework, the conventional method to compute it involves extracting the transverse-traceless (TT) parts of the tensor perturbations that are sourced by the matter fields  \nρ GW = 1 32 πG ⟨ ˙ h TT ij ˙ h ij TT ⟩ , h ij +3 H ˙ h ij -1 a 2 ∇ 2 h ij = 16 πGS ij , (44)  \nand can be projected as  \nh TT ij ( t, x ) = ∫ d 3 k (2 π ) 3 e i k · x Λ ij,kl ( ˆ k ) h lm ( t , ˆ k ) , (45)  \nwhere the projector operator is defined as Λ ij,kl ( ˆ k ) ≡ P il ( ˆ k ) P jm ( ˆ k ) -P ij ( ˆ k ) P lm ( ˆ k ) / 2, with ˆ k = k /k and P ij ( ˆ k ) = δ ij -ˆ k i ˆ k j .  \nIn NR simulations, GWs are often extracted by projecting ˜ A ij into the TT gauge  \nρ GW = 1 8 πG ⟨ ˜ A TT ij ˜ A ij TT ⟩ , (46)  \nNote, however, that although the GW spectra can be directly extracted by projecting out the transversetraceless part of the extrinsic curvature, this is only approximate as it assumes ˙ h ij ≈ -2 ˜ A ij , which may not be true when perturbations are large. The GW spectrum Ω GW ≡ ρ -1 c d ρ GW / dlog k can be computed from lattice simulations of length L as [102]  \nΩ GW ( k ) ≈ 3 2 π 2 L 3 k 3 ⟨ K ⟩ 2 ∑ i,j | ˜ A ij ( k, τ ) | 2 . (47)  \nKou et al. [98] followed a similar approach defining an expanding background metric to be subtracted from the full metric accounting for the conformal factor ψ , with δγ ij = γ ij -⟨ ψ 4 ⟩ δ ij . The energy density in GWs was then computed as  \nρ GW = 1 32 πG ⟨ δγ TT ij ;0 δγ ij TT ;0 ⟩ , (48)  \nwhere ; 0 is the covariant derivative compatible with the background metric. In Fig. 7 we show their results, where they compared the extracted Ω GW from oscillon preheating using a rigid FLRW scheme and NR simulations with different gauges: (i) synchronous gauge (geodesic slicing): α = 1 and β i = 0; (ii) A combination of the 1+log slicing and Gamma driver shift conditions, Eqn. (40); (iii) Radiation gauge in Eqn. (41). This comparative approach aimed to understand the impact of NR gauge choices on the extraction of gravitational waves in a fully relativistic setting, and they concluded that the extracted GWs are highly sensitive to the choice of gauge. Synchronous gauge (geodesic slicing) was found to be an unreliable gauge condition for extracting GW spectra, as it contained redundant gauge modes that artificially enhanced the GW power spectrum by 1-2 orders of magnitude.  \nThis example highlights the importance of carefully choosing and interpreting diagnostics from NR simulations.  \n<!-- image -->  \n<!-- image -->  \nFig. 7 : Time evolution (from left to right) of GW spectra for a preheating model with different NR gauge choices. The figure shows the breakdown of the synchronous gauge in reliably extracting the GW spectra. The insets in the bottom panels depict the ratio between the GW spectrum in the corresponding gauge condition and that of the FLRW simulation. Figure from Kou et al. [98].  \n<!-- image -->", '2.7 Units': 'Cosmologists use units in which ℏ = c = 1, and quantities are described in terms of the reduced Planck mass M Pl = √ ℏ c/ 8 πG . In this system of units masses ∼ M Pl and lengths ∼ 1 /M Pl , unlike in NR where G = c = 1 and both lengths and masses have units of [ M ] where M is some physical mass scale in the problem (e.g. M = M ⊙ ). So, for example, energy density has units of [ M Pl ] 4 in cosmology and [ M ] 2 in NR. To convert a quantity X with units of [ M ] n in NR units to its cosmological value Y in units of [ M Pl ] m , one uses  \nY = X ( M/M Pl ) n . (49)  \nIn most simulations we are dealing with length scales well above the Planck length 10 (where we expect the classical description to hold), and so ( M/M Pl ) ≫ 1.  \nIn NR, there remains a free mass/length scale M in the problem that once set determines the physical values of the results. This is usually set such that the code values are order 1 numbers, which is advantageous for numerical reasons. For example, in an equal-mass binary BH simulation, the BH bare masses are chosen to be M ∼ 0 . 5, with exact values set such that the total ADM mass is exactly 1. This simulation describes equally well a merger of supermassive BHs and solar mass ones, it simply requires the results to be scaled appropriately between the two. So in NR papers one sees quantities like lengths and times described in units of M , which in this example is the total ADM mass. The fact that the BH mass and not its Schwarzschild radius was chosen as the NR reference scale is an unfortunate historical convention that causes a lot of unnecessary confusion for those coming from Planck units. It is much better to think of the fundamental scale as being set not by the BH mass, but rather by the Schwarzschild radius relating to the curvature created by such a mass. Consider for example the mass parameter µ that appears in the Klein-Gordon equation  \ng µν ∇ µ ∇ ν ϕ = µ 2 ϕ . (50)  \nIn an NR simulation of a scalar field and a BH, the physical scale of interest in this system could be set by either µ , or the mass M BH of the BH; once one is chosen, the other is fixed. Say that we choose to set the code values M BH = µ = 1. This does not mean that if we choose the BH to represent a mass of M ⊙ then the scalar particle mass is M ⊙ , which is clearly nonsense. It means that the Compton wavelength of the particle is of the same order as the Schwarzschild radius of the BH, which is about 1km. The correct particle mass m p can be obtained by recognising that µ = m p c/ ℏ in geometric units does not have units of [ M ] but rather [ M ] -1 (or more sensibly, inverse length), so applying the formula in Eqn. (49) (with m = 1 and n = -1)  \nµ Planck = µ Geo ( M ⊙ M Pl ) -1 M Pl = 1 ( 2 × 10 30 kg 2 × 10 -8 kg ) -1 × (10 28 eV)  \n̸  \nwhich one can verify corresponds to particle with a Compton wavelength of order ∼ km. The key takeaway here should be that in NR simulations ℏ = 1 and this not only changes the values but also the dimension of the quantities. It is important to go back to the real physical unit (in terms of mass, length and time) and check the dimension in each unit system to get the right transformation, with [ M ] ∼ [ L ] ∼ [ T ] in geometric units and [ M ] ∼ 1 / [ L ] ∼ 1 / [ T ] in Planck units.', '3 Genesis: The early universe': "We start our review of cosmological NR works with those studying the beginning of our universe, or (depending on one's preferred paradigm) those describing previous epochs before the standard hot Big Bang cosmology commences. The broad goal of these works is to understand the nature of the cosmological singularity, to connect it to our currently observed homogeneous and isotropic universe, or to remove it entirely by postulating a bounce in the cosmological evolution that would see the universe transition between collapsing and expanding phases.  \nMost of these works have a common theme in that they seek to understand how generic and robust certain processes are within the context of classical GR in cosmological spacetimes. We can only observe one Universe, and so in this context NR provides the potential for numerical experiments to understand what other possibile universes may have been generated if we assume a certain paradigm. Such studies are open to philosophical arguments about anthropics or how to precisely define terms like generic or robust (in particular, the so-called measure problem , see for example [103]). Nevertheless, they can at the very least provide material with which to inform the debate, and at most they provide strong statements about consistency in models invoked to explain the observed features of our universe. In the final subsection we review work to identify potential signatures from inflationary bubble collisions, which provides a good example of attempts to connect such early universe paradigms with concrete observables.", '3.1 Cosmological singularities and mixmaster universes': "Some of the earliest work in NR was applied to understanding the nature of singularities (see [6] for a review). Here we focus on aspects related to the cosmological singularity case, which has also been discussed in the reviews [5-8, 10-12]. Note that here when we discuss time evolution, usually the evolution is running backwards towards the singularity from some initial state that represents our universe at a later time. As a result the spacetimes are collapsing and terms that normally dilute with expansion may become dominant, such as curvature and anisotropy.  \nThe fact that rewinding our universe unavoidably results in a singularity in the absence of exotic matter types is a consequence of Penrose's singularity theorems [104]. However, these theorems simply tell us that null geodesics cannot be extended beyond some future event, and not about the nature of the singularity that is reached. Belinskii, Khalatnikov, and Lifshitz (BKL) conjectured that the singularity in generic gravitational collapse is space-like, local, and oscillatory [105]; more specifically, that at each point they approach a state that is well-described by the vacuum, homogeneous Bianchi type IX (Misner's Mixmaster [106]) cosmology.  \nThe statement that the singularity is space-like means that it cannot be observed before it is reached. The local aspect of the spacetime is crucial in ensuring that the dynamics of neighbouring points become decoupled (i.e. spatial gradients become negligible in terms of their physical effects on the local evolution). In the resulting Mixmaster cosmology, the spacetime transitions through a series of Kasner spacetimes for which the metric is  \ng µν dx µ dx ν = -dt 2 + t 2 p 1 dx 2 + t 2 p 2 dy 2 + t 2 p 3 dz 2 (52)  \nwhere the exponents p i are constants (a spatially flat FLRW cosmology is a special case where the values of p i are equal). In the vacuum case the exponents must satisfy ∑ i p i = ∑ i p 2 i = 1, which implies that one direction must expand whilst the other two contract. The transitions between epochs with different Kasner exponents are discrete and chaotic, but have a well-defined map, see Barrow [107]. They occur when the local spatial curvature (3) R becomes dominant (as it may in a closed FLRW cosmology or where the anisotropy of the spatial metric grows during collapse). In cosmologies where the (3) R never dominates the spacetime will remain stuck in a particular Kasner case. This is then referred to as an Asymptotically Velocity Term Dominated (AVTD) singularity since only the kinetic term in the Hamiltonian drives the  \n<!-- image -->  \nFig. 8 : (Left:) Spike development in the metric quantity N α α (the quantity that vanishes on a spike surface) in a 1+1D simulation. (Right:) Spikes in N α α in a 2+1D simulation manifest as curves. From Garfinkle & Pretorius [125].  \n<!-- image -->  \nevolution at late times. The BKL conjecture is that the Mixmaster dynamics are the most generic case, with AVTD spacetimes only occurring where particular symmetries are imposed. Therefore the nature of the generic singularity is of an oscillatory Mixmaster type.  \nEvidence in support of the BKL conjecture has been obtained from several NR works. The first studies by Berger & Moncrief [108] for T 3 × R vacuum Gowdy spacetimes that have two spatial symmetry directions found only AVTD behaviour, as did later extensions to Gowdy spacetimes with one spacelike U(1) symmetry or other topologies [109-112]. However, works with more generality, and including matter, were found to be consistent with a BKL picture [113-118]. The general vacuum case with T 3 × R topology but no spatial symmetries was first studied by Garfinkle in [114] using a scaled tetrad formalism of Uggla et al. [40]. The connection of BKL singularities to the null singularities expected in generic BHs or other asymptotically flat spacetimes is still under investigation [119, 120].  \nWhilst numerical evidence has mostly supported the BKL conjecture, the initial work by Berger & Moncrief [108] was also responsible for the discovery of the existence of a new feature of small-scale structure, now called spikes . These occur where points are stuck in an old Kasner epoch while points around them transition to new ones. The term spikes is a bit of a misnomer because these are generically 2D surfaces in a 3+1D spacetime, but since they are codimension one surfaces in a spacetime with symmetry, and were mostly studied in 3+1D spacetimes with 2 symmetry directions, they then appear to be point-like in plots like the left panel of Fig. 8. Studies of spikes have been performed in G 2 cosmologies [121] and Gowdy spacetimes [122-124]. Recently adaptive mesh refinement (AMR) has been used by Garfinkle & Pretorius to permit the study of cases with only one symmetry direction (so that the spikes manifest as closed curves [125], see Fig. 8. The numerical simulations are found to be consistent with an analytical formula for the shape of the spikes that suggests that whilst the gradients at spikes may grow arbitrarily large, their impact on the dynamics nevertheless decreases as the singularity is approached. This then implies that the spikes may be only transient features, and that the BKL picture holds even at these somewhat special points.  \nThe BKL conjecture is also sometimes extended with the statement that matter does not matter , implying that the vacuum dynamics are sufficient to describe the approach to the singularity. This is evident for most regular matter types, for which the the growth in the energy density during collapse is subdominant to the growth in anisotropy (which grows as ∝ 1 /a 6 ). However, one exception is a scalar field in which the energy density may grow at the same rate, or greater depending on its potential, which can suppress the BKL behaviour [39, 126-128]. This fact is used as a smoothing mechanism in bouncing scenarios, as will be discussed in Sec. 3.3.  \nThe simulations used to study singularities inspired a number of novel numerical techniques, for example the tetrad formulation that permits the overall scaling of the expansion to be factored out [40, 41], tracking swarms of test particles to measure gravitational waves [129] and symplectic integration schemes for accuracy and stability [130, 131]. The use of GHC by Garfinkle [39] inspired the successful  \nsolution to the binary BH problem by Pretorius [37, 38], nicely demonstrating how results from different areas of NR can find applications in other physical scenarios.", '3.2 Inflation and the initial condition problem': 'Cosmic inflation [132-136] is a theory of the early universe in which a period of nearly exponential expansion occurs, providing an explanation of why the universe on large scales is spatially flat and homogeneous, despite the expectation that the cosmological singularity will be chaotic and, due to its local nature, uncorrelated at separated points. Inflation also provides a natural mechanism to generate an approximately scale-invariant spectrum of primordial perturbations that seed structure formation, as observed in the CMB [137]. The simplest mechanism for inflation is a single scalar field, the inflaton, that slowly rolls down a roughly flat region of its potential. To be consistent with this picture, sufficiently generic (presumably inhomogeneous) initial data should still permit the scalar field potential to dominate and at least some of the universe to undergo a period of inflation. Single-field models of inflation are classified as high-scale or low-scale models, depending on the energy scale associated with inflation, or as large-field or small-field models depending on the field range over which the minimum required 60 e-folds of inflation occurs (some care is needed in interpreting the different terms [138]). In cases where the inhomogeneities are strongly backreacting at the Hubble scale (which is predominantly the case for large-field models), NR simulations may be needed to determine whether this happens.  \nThe problem of initial conditions for inflation has been studied using many analytic and semi-analytic methods (see Brandenberger [139] for a review). Whilst anisotropy was shown analytically to remain subdominant in work by Turner & Widrow [140], the general case of strong inhomogeneities required a numerical treatment. Using NR, early work by Goldwirth & Piran [81, 141, 142], focused on spherically symmetric inhomogeneities in an inflationary theory, assuming a closed universe to avoid the issue of imposing boundary conditions. The initial work included a massless scalar or fluid radiation component along with the inflaton scalar field, the latter of which was subject to a range of potentials. These simulations suggested that inflation would begin as long as there is initially a homogeneous region of size R ∼ H -1 inf at the inflationary part of the potential, where H inf is the inflationary scale.  \nFor high-scale inflation, a number of simulations have supported the claim that inflation persists with generic initial conditions. The first 3+1D simulations with periodic boundary conditions by Laguna, Kurki-Suonio & Matzner [143, 144] showed that inflation could occur even if kinetic and gradient energies were initially comparable to the potential energy density. However, these simulations were limited to initial conditions that did not lead to BH formation. Simulations leading to BH formation were carried out by Shibata et al. [54] for axisymmetric spacetimes with a cosmological constant and an initially large distribution of gravitational waves. While their boundary conditions (asymptotically de Sitter) and the presence of a cosmological constant meant that the end state is always de Sitter [145], they showed that there can be an initial gravitational wave radiation dominated period (i.e. the GW do not always collapse to form a Schwarzschild-de Sitter spacetime.) It would therefore be interesting to revisit the work with a more general inflationary model and/or boundary setup. More general 3+1D initial conditions, including those that lead to BH formation, were considered first in East et al. [90] and studied later by others [72, 75, 80, 91-93, 146]. These showed that high-scale inflation can occur in the presence of large inhomogeneities even where BHs form, see Fig. 9. Some recent work by Garfinkle, Ijjas & Steinhardt has opposed this conclusion on the basis that the initial conditions and boundary choices made may bias the results towards successful inflation in some non-generic regions of low initial curvature [49, 147], see Joana [148] for a partial response.  \nThe robustness of low-scale models to inhomogeneous initial conditions is an even more important issue. Since the energy scale is smaller, the size of the causally disconnected patches is larger, and the potential energy may be orders of magnitudes smaller than other components. Consequently, it is possible that large inhomogeneities could spoil inflation before it begins. This has motivated a wealth of studies on the effects of initial conditions on low-scale inflation [72, 75, 90-93]. In general, low-scale models do seem to be less robust to inhomogeneities, with the field dynamics often disrupting the inflationary slow-roll. An important observation by Aurrekoetxea et al. [93] is that the convexity of the inflationary potential plays a strong role in its robustness, with convex models naturally more robust. This can be understood from arguments about the tendency of the field gradients to homogenise the field, in competition with the tendency of the potential gradient to enhance them.  \nThe effect of large non-uniform initial inflaton velocity in addition to an inhomogeneous initial field profile was first studied in Corman et al. [72] for the large-field case in both single and two-field inflationary models, and it was found that these models are generally robust, a conclusion that was later  \nFig. 10 : (Left:) Figure from [72] showing gradient (blue solid), kinetic (red dashed) and potential (green dot-dash) contributions to the average energy density. Despite initially high kinetic and gradient contributions, the potential dominates at late times, permitting inflation to occur. Gray regions indicate black hole formation. (Right:) Figure from [75] showing the spread in the resulting number of e-folds in the case of a single mode of initial field perturbations, with homogeneous (blue), single mode (red) and double mode (green) perturbations in the kinetic term. A homogeneous kinetic term has the most significant impact, and can significantly reduce the number of e-folds of inflation. In other cases some parts of the spacetime still always experience significant inflation.  \n<!-- image -->  \n<!-- image -->  \ntHo  \nFig. 9 : (Left:) Figure from East et al. [90] showing the first successful simulations of high scale inflation in the 3+1D case in which BHs form. (Right:) Figure from Clough et al. [91] showing the mass of the BHs formed versus the size of the initial perturbations in a large field model, showing that the maximum mass of the BH that can be formed by such perturbations is bounded.  \n<!-- image -->  \n<!-- image -->  \n〈N〉  \nsupported by the work of Elley et al. [75] using different methods for generating the initial conditions. In Corman et al. [72], the XCTS method with a spatially constant mean curvature K and a rescaled scalar field ansatz were used to generate initial data, as well as an ansatz for the momentum constraint suggested by Garfinkle & Mead [67], whereas the latter paper used the CTTK method of Aurrekoetxea et al. [48] that creates a spatially varying mean curvature K and greater gradients in the extrinsic curvature A ij . Despite the different initial data choices, and the use of GHC versus BSSN formulations, both works showed that the effect of velocity on low-scale models resulted in a larger spread of e-folds in different regions [75], but did not significantly improve or harm the success of such models, in comparison to the zero velocity case - see Fig. 10. This work addressed the concerns of Ijjas [13] regarding the effect of kinetic inhomogeneities (but does not respond to the more recent criticisms in [49, 147]).  \nBeyond Corman et al. [72], little work so far has focussed on the multi-field case. Joana [149] investigated the robustness of the Higgs inflationary model to inhomogeneous multi-field initial conditions,  \n˜  \nη  \nwith a non-minimal coupling to gravity, as well as the subsequent preheating period, showing that the dynamics from gravitational shear and tensor modes can only delay, but not prevent, cosmic inflation. Given that higher energy justifications for inflationary potentials (e.g. string theory) would generically predict a multi-field landscape, more work on the robustness of such models is justified.  \nThe impact of high-energy modifications to GR are also an interesting direction of study, given that at the start of inflation such corrections may well have been relevant. Recent advances in well-posed formulations of more general 4-derivative models like Einstein-scalar-Gauss-Bonnet [150] and higher derivative pure tensor models [151] have opened up the possibility of studying the impact of higher order EFT corrections in such models. Work to develop the initial data methods, and evolutions with a single field acting both as the inflaton and the additional scalar degree of freedom have been completed by Brady et al. [76, 152]. Assuming that inflation happens at all with such models, the additional effect of inhomogeneities via the derivative terms is not found to be significant. In particular, the dynamics remain in the weakly coupled regime, rather than being driven out of the EFT regime, and the dominant effect is the additional tilt of the effective potential created by the change to the background/average cosmology. Models with a separate inflaton field and scalar degree of freedom, and pure tensor higher derivative models may give rise to more interesting dynamics and will be a focus of future work.  \nFinally we note that most of the work described here is quite binary in its conclusions - inflation either works or it does not. Braden et al. [153] provides an interesting example of using NR simulations to provide actual observational constraints, assessing the CMB quadrupole constraint on the amplitude of the initial perturbations and the size of the observable universe relative to a length scale characterising ultra large-scale structure. Such studies are challenging because of the cost of NR simulations, but this study used a toy model reduced to 1D perturbations to make progress.', '3.3 Ekpyrotic and bouncing scenarios': 'The ekypyrosis model of cosmology [154] posits that the hot Big Bang is generated by the collisions of a brane and its bounding orbitfold in the bulk (extra-dimension) space. During the infall phase of the brane collision, the 3D universe on the brane (on which we live) is in a matter-dominated phase and undergoes a slow contraction. This phase ends at the collision, and the universe reheats during the collision and bounces, reentering the hot Big Bang phase as it falls away from the orbitfold. This phase cools, until the brane turns around and undergoes infall again, resulting in a cyclic universe with repeated bounces.  \nThis string-inspired model has developed into a more general paradigm in which the universe undergoes cyclic periods of expansion and collapse, separated by a non-singular bounce [155, 156]. In such cases the slow contraction phase leading to smoothing of inhomogeneities, rather than an exponential expansion as in inflation (Ijjas & Steinhardt provide a pedagogical review in [157]).  \nIt was suggested that the slow contraction phase could be an efficient smoother , i.e. it could drive the universe from an inhomogeneous phase into a homogeneous one [158]. The underlying mechanism for this smoothing is ultralocality [159], in which during the collapse the universe evolves into a phase where the gradients become much less important than the kinetic term, and hence the set of PDEs become a set of ODEs as in the BKL conjecture of [105], discussed in Sec. 3.1. To generate such an evolution, one can consider scalar matter with a negative-definite potential V ( ϕ ) < 0. In a homogeneous FLRW background, the equation of state parameter w is then  \nw = (1 / 2) ˙ ϕ 2 -V (1 / 2) ˙ ϕ 2 + V . (53)  \nand one can see that if V < 0, the in general w ≥ 1 as long as | V | < (1 / 2) ˙ ϕ 2 . Meanwhile, the Hubble parameter for an FLRW spacetime is given by  \n| H -1 | ∝ a ϵ , (54)  \nwhere ϵ = (3 / 2)(1 + w ). In standard cosmological evolutions ϵ = 2 and ϵ = 3 / 2 for radiation and cold dark matter domination respectively, inhomogeneities grow. However, in slow contraction ( ϵ > 3), the scale factor shrinks much more slowly than the Hubble radius, and hence the universe gets driven to an ultralocal phase which homogenises it 11 . Clearly, this argument assumes a homogeneous universe to begin with, and hence numerical methods must be employed to test this hypothesis.  \n<!-- image -->  \nFig. 11 : (Left:) Figure from the first work on slow contraction by Garfinkle et al. [160]. The time evolution of the normalized energy densities in matter (solid), curvature (dot-dash) and shear (dashed), showing that some regions become homogeneous and isotropic. These regions eventually dominate the proper volume. (Right:) Figure from Cook et al. [161], recreating the original simulations ((a) top row) and showing that the inhomogeneous regions disappear if one chooses the scalar to be rolling down the potential at all points initially ((b) bottom row), which the authors argue is more physically motivated.  \n<!-- image -->  \nThe first NR work to simulate slow contraction by Garfinkle et al. [160] employed the tetrad formalism of Uggla et al. [40, 162]. Using a scalar potential of the form  \nV = -V 0 exp( -kϕ ) , (55)  \nthey simulated the evolution of slow contraction spacetimes with initial inhomogeneities in the scalar and extrinsic curvature terms in a CMC gauge (see Sec. 2.5) with perturbations in a single spatial direction. They showed that, for some initial inhomogeneous configurations at least, the domain bifurcate into homogeneous w ≫ 1 slow-contraction regions and an inhomogeneous w ≈ 1 regions, with the ratio of proper volume of the former to the latter regions growing exponentially fast. In particular, in the slow-contraction region, the space becomes smoother, as expected. See Fig. 11 left panel.  \nThis work was later extended by Cook, Ijjas and others [161, 163], who argue that the initial condition for the scalar field in Garfinkle et al. [160] was inconsistent with starting slow contraction. Specifically, the scalar field must roll towards increasingly negative V in Eqn. (55), and if the initial conditions are set such that this is true initially at all points in the initial hyperslice, then all the regions generically enter the slow-contraction phase that exhibits efficient smoothing, see Fig. 11 right panel. In Kist & Ijjas [164], it was shown that the slow contraction phase is generic and does not depend strongly on the shape of the potential, as long as it is sufficiently steep M Pl | V, ϕ/V | ≥ 5, while in Ijjas et al. [165] it was shown that these results are robust to perturbations in two spatial dimensions and additional perturbations modes.  \nIn addition to smoothing, one of the key requirements for any model of cosmogenesis is the generation of a scale invariant spectrum of perturbations. Unlike inflationary models, in which the perturbations of the inflaton can provide this spectrum, slow contraction models require another scalar field χ which sources the scale-invariant perturbations, with the following non-minimal coupling to the kinetic term of the ϕ field [166, 167]  \nL int = -1 2 e -ϕ/m ∇ µ χ ∇ µ χ , (56)  \nfor some mass scale m ≤ M , where M is the mass of ϕ . In Ijjas et al. [168], this model was numerically studied and it was found that, unlike in the single field case, the presence of the χ field can destabilise the system as ϕ energy can leak into χ . Nevertheless, they demonstrated that this timescale of instability is longer that the bounce cycle time as long as the mass of χ is small but non-zero.  \nFinally, while the above works focus on the smoothing physics of the slow contraction phase of a bouncing cosmology scenario, Xue et al. [169] focused on the physics of the bounce itself. To execute the bounce after a slow contraction phase, an additional massless scalar ghost field χ (i.e. a scalar field with wrong sign kinetic term) was added  \nL = -1 2 ( ∂ϕ ) 2 -V ( ϕ ) + 1 2 ( ∂χ ) 2 , (57)  \n<!-- image -->  \nFig. 12 : (Left:) Figure from Ijjas et al. [168] showing slow contraction in the multi-field case needed to generate scale invariant perturbations. Time evolution of the normalized energy densities in the fields (green) the curvature (red) and the shear (blue) are shown. The model drives the universe towards FRW, but a non-zero gradient in the auxilliary field destabilizes the FRW state on longer timescales. (Right:) Figure from Xue et al. [169] showing the local expansion as a function of space, evolving over time. Whilst the whole of the space bounces (going from negative to positive expansion), the space in the middle bounces at a much later time compared to the other regions, creating some inhomogeneity at late times.  \n<!-- image -->  \nto the usual slow contraction field ϕ . This ghost field violates the Null Energy Condition, but since its energy density grows negative faster than the positive energy of ϕ as the scale factor a ( t ) → 0, it will dominate at small a ( t ) and turns the contraction phase into an expansion phase, executing the bounce [170]. Xue et al. [169] studied the evolution of perturbations through the bounce using the harmonic gauge (see Sec. 2.5), and found that if the initial perturbations are large, the induced non-linearities can cause inhomogeneous regions to collapse into a singularity while relatively homogeneous regions will undergo the bounce successfully. In a slightly different direction Corman et al. [53] studied the behaviour of BHs through such a two-field bounce model, and found that for sufficiently large BHs (with horizons ∼ 3 -4 times the order of the minimum cosmological horizon) the two horizons can merge and disappear for a short time during the contracting phase. Nevertheless, this state of affairs is temporary and as the bounce enters the expanding phase, the horizons separate and the BH end state is similar to its initial pre-contraction phase value.', '3.4 Signatures of inflationary bubble collisions': 'When a spacetime region residing in a metastable vacuum transitions quantum mechanically to a lower energy state, via the nucleation and expansion of a bubble of the new phase [171, 172], eternal inflation can occur [173-175]. In such a case our observable universe is just one of many bubble universes within an infinite inflating space. If the rate of nucleation is lower than the expansion rate of the parent vacuum, the original phase is never consumed, leading to eternal inflation (see Refs. [176, 177] for a review). Collisions between expanding bubbles could leave distinct observational signatures in the CMB [178, 179], which may allow us to confirm this hypothesis with data. Due to the non-linear amd dynamical nature of the problem, numerical simulations are needed.  \nThe spacetime of two colliding bubbles possesses SO(2,1) symmetry, arising from the intersection of the hyperbolic SO(3,1) symmetry of individual bubbles. This allows the collision to be treated as a 1+1D problem. Early numerical simulations have studied the dynamics of isolated vacuum bubbles [180] and collisions between them [181-186] (particularly in the context of branes). However, these bubbles lack an interior cosmology and hence are not adequate for making phenomenological predictions. Johnson et al. [187] performed the first extensive study of the outcome of cosmic bubble collisions using 1+1 dimensional NR simulations. They focused on the collision between thin-wall bubbles for potentials with different features and compared their simulations to several commonly used approximations. Where the interior of the expanding bubble was vacuum, their results showed good agreement with the Israel Junction Condition formalism [188] that matches GR solutions across a boundary, but this did not hold when the bubbles contained an interior cosmology. In that case, they confirmed the validity of the free-passage approximation [189-192], in which the field profiles of the colliding bubbles are superposed. The authors concluded that the structure of the potential was more important for determining the outcome of a collision than the kinematics, with most of the energy going into the formation of oscillons and pockets  \nFig. 13 : Schematic figure depicting the procedure from a scalar inflationary potential to a prediction of the observable comoving curvature perturbation R . NR simulations are performed to construct the full collision spacetime (upper right). The perturbed FRW metric inside the observation bubble is reconstructed by evolving a set of geodesics through the simulation (lower right panel). A gauge transformation is used to extract the comoving curvature perturbation R late in the inflationary epoch. Evolving the comoving curvature perturbation, the observable universe is split into a region that is, and a region that is not, affected by the collision. The lower left panel depicts the surface of last scattering inside the observation bubble in a reference frame where the observer is at the origin of coordinates. Figure from Wainwright et al. [194].  \n<!-- image -->', 'Cosmological Signatures': 'of false vacuum, as opposed to scalar radiation. They found that bubble collisions could disrupt inflation if it was driven by a small-field model, but were less likely to do so for large-field inflationary models.  \nIn a series of papers, Wainwright et al. [193-195] made the first quantitative connection between the scalar field inflationary model and the signature imprinted on the CMB - an azimuthally symmetric temperature anisotropy. The anisotropy can be computed in the Sachs-Wolfe limit as  \nδT T = R ( ⃗ X ls ) 5 , (58)  \nwhere R ( ⃗ X ls ) is the comoving curvature perturbation on the surface of last-scattering. In the first of the papers [193], the authors studied the setup depicted schematically in Fig. 13 to obtain R in the observer frame. They found that the comoving curvature perturbation is well-described by a power law, with its shape and amplitude determined by four phenomenological parameters: the total number of predicted collisions, the shape of the power law, the radius of our past lightcone at last-scattering, and the angular size.  \nThe signatures depend on several features that depend on the scalar potential: the post-nucleation slow-roll inflationary period, the post-inflationary cosmology, the properties of the cosmology in the colliding bubble, and the details of the potential barriers surrounding the false vacuum. In the second paper of the series [194] the authors focused on models that predict only one observable collision and studied the effect of varying the potential barrier in collisions between identical and non-identical bubbles. The collision between identical bubbles led to an advance of the inflaton, a positive curvature perturbation R > 0 and thus a hot spot in the CMB. The collision between non-identical bubbles, which tunnel in opposite directions, caused a delay of the inflaton, R < 0, and resulted in a cold spot instead.  \nThe earlier works were limited to extracting the signatures in the vicinity of the collision boundary, and could not compute the observables far inside the collision region. In the final paper [195], the authors extended their techniques to explore the full bubble-collision spacetime, allowing them to make predictions for all observers. In the case of identical bubbles they found that observers in the instantonand overlap-born region near the collision boundary see different signatures, and oscillatory contributions to the curvature perturbation for non-identical bubbles.  \nA similar numerical framework was used by Johnson & Lin [196] to compute the observational signatures of a classical transition in which a collision between bubbles produced a region with a new vacuum, different from either of the vacua in the colliding bubbles. They found that the main observables were negative spatial curvature and an approximately quadratic comoving curvature perturbation with planar symmetry, which again could be mapped to a temperature quadrupole in the CMB.  \nAn alternative approach to studying bubble collisions with gravitation beyond the thin-wall approximation is the use of the double-null formalism [197, 198], which consists in decomposing the spacetime interval in the null ( u, v ) coordinates. Hwang et al. [199] used these techniques to study the collision between vacuum bubbles with spherical, planar and hyperbolic symmetry. Kim et al. [200] followed a similar approach to compute the gravitational waves emitted from the collision using the quadrupole approximation (see Sec. 2.6). They found that the field value of the false vacuum controlled both the frequency of the gravitational waves and the amplitude of the waveforms. However, due to the highly relativistic nature of the bubbles when they collided, the results obtained using the quadrupole approximation should be considered qualitative rather than precise predictions.', '4 Exodus: The transitional universe': 'The era of Big Bang Nucleosynthesis (BBN) at around T = 1 MeV marks the beginning of the period in which we can say that (broadly speaking) we understand how the universe has evolved on large scales. Before this era, it is here-be-dragons observation-wise until the epoch in which the primordial spectrum of perturbations are seeded. The fact that the universe is expanding, and that BBN is hot, means that at some point in the far past the universe must either begin hot, or have undergone an epoch of reheating after the perturbations were seeded. In particular, in the inflationary model of cosmogenesis, the reheating process can feature non-perturbative dynamics that give rise to exotic compact objects and a background of GWs. Given that the early universe was opaque to light before the release of the CMB, such primordial GWs may be the only way to directly probe such early epochs, so characterizing their amplitude and spectral shape is important.  \nAs the universe cooled, it is expected that it underwent phase transitions that broke down a unified physical model into the forces that govern our universe today. During this process various topological defects may have arisen, such as domain walls, cosmic strings, and monopoles, depending on the symmetry of the vacuum. These objects can leave imprints on the CMB and in GW backgrounds. Primordial black holes (PBHs) may also have formed during the early phases of the universe from the collapse of large density fluctuations generated during inflation or as the end products of phase transitions and dynamics of topological defects. In some mass bands, PBHs could constitute all (or a substantial part) of the dark matter of the universe and provide a mechanism for the origin of supermassive BHs.  \nMany simulations of the early universe neglect the local effects of gravity, accounting only for the overall FLRW expansion. Whilst in many cases such an approximation is justified, in cases where high density contrasts are reached, in particular in PBH formation, NR simulations will be needed. In this section, we will review the results from cosmological NR studies of reheating, topological defects and PBH formation.', '4.1 Reheating and oscillon formation': 'After inflation, the universe is dominated by the scalar field that drives the exponential expansion - the inflaton. At some point, the universe must transition between the inflaton dominated era to a radiation dominated universe, converting the energy density in a process known as reheating . This mechanism sets the scene for the hot Big Bang, converting the vacuum energy in the inflaton into radiation via its decay into other lighter fields. However, a universe filled with an oscillating inflaton condensate is governed by Mathieu-like equations, which can, depending on the inflationary potential, grow the field amplitudes exponentially and drive explosive particle production via parametric resonance. This non-perturbative process is known as preheating (see Amin et al. [192] for a review).  \nFig. 14 : (Left:) The power spectrum of the gauge-invariant metric perturbation Φ is plotted against k , with all modes initially excited. Each panel corresponds to a different time. (Right:) The growth of C 2 /R 2 during the simulation. Figures from Easther & Parry [82].  \n<!-- image -->  \nPreheating takes place far from thermal equilibrium and is an inherently non-linear process. Several levels of approximation can be used: (1) Linear perturbation theory: The field amplitudes are assumed to be small so that mode-mode couplings can be ignored. (2) Lattice field theory on a rigid FLRW background: Mode-mode couplings are allowed but the backreaction onto spacetime is neglected. (3) Lattice field theory with GR: The full dynamics of the field and the spacetime are evolved. The final case can be used to quantify the impact of the assumptions made on the conclusions reached in the first two cases. Some of the earliest preheating studies that included GR focused on comparing these three cases. Gravitational effects from larger inhomogeneities may lead to collapse on smaller scales, and on scales much larger than the Hubble length, non-linear effects might enhance metric perturbations, with potentially observable impacts on the CMB.  \nParry & Easther [201] pioneered numerically solving Einstein field equations in 1+1D with planar symmetry and periodic boundary conditions to study preheating after the end of m 2 ϕ 2 inflation. The simulations showed that wavelengths considerably larger than the Hubble length k ≪ < 1 /H would not undergo parametric amplification for the m 2 ϕ 2 model, thus verifying the perturbative treatment of the post-inflationary era in the absence of a resonance [202]. However, by studying perturbations with larger k they demonstrated that non-linear effects can transfer power between resonant and non-resonant modes. In over-dense regions the onset of gravitational collapse was observed and they concluded that non-linear gravitational effects would broaden any resonance band seen in perturbative analyses. This was investigated further by the same authors for the λϕ 4 model [82]. In the case of a single excited mode, the perturbative evolution initially agreed with the non-linear treatment, However, while perturbation theory predicted indefinite exponential growth, non-perturbative effects in the field dynamics eventually terminated the resonance. Comparing the rigid background and non-linear gravity simulations showed that when self-gravity was included, the amplitude of perturbations remained comparatively large, whereas without it they decayed. When multiple inhomogeneous modes were initially excited, power was transferred to long wavelength modes in both the rigid FLRW background and the full GR simulations, see Fig. 14. In the GR case, the overall growth was found to be larger and began sooner, although the authors could not conclude whether this was due to backreaction effects or other choices made for the simulations, such as periodic boundary conditions. To study the impact of gravity, they suggested tracking the ratio of the domain integrals of the Weyl curvature scalar C 2 ≡ C µνγλ C µνγλ and the usual Ricci curvature scalar, C 2 / R 2 , as a measure of how far the spacetime departs from perfectly FLRW, see Fig. 14. As they did not observe this quantity diverging, they concluded that no PBH formation occurred for the λϕ 4 theory. This is in contrast to Finelly & Khlebnikov [203], who performed general relativistic 1+1D simulations of preheating following hybrid inflation in spherical symmetry. They found that in this scenario, super-Hubble metric perturbations could experience significant growth, with potential evidence for BH formation.  \nThe first 3+1D studies of preheating by Bastero-Gil et al. [83, 204, 205] employed the ADM formalism with the synchronous gauge to study m 2 ϕ 2 + λϕ 4 and hybrid models of inflation. For the chaotic model they found that the amplification of superhorizon modes observed in 1+1D simulations [82, 203] persisted in 3+1D. For the hybrid model, where the resonance is more explosive, scalar metric perturbations were found to have little impact on the evolution of field fluctuations. They compared the obtained energy  \nFig. 15 : Oscillon formation during preheating. The faces of the box show the evolution of inflationary fluctuations and the inner regions correspond to the isocontours of density, illustrating the formation of oscillons. Figure produced with simulations from Aurrekoetxea et al. [73].  \n<!-- image -->  \ndensity in gravitational waves to integrating the tensor modes with a rigid FLRW scheme and found that it was approximately an order of magnitude larger 12 .  \nMore recently, Giblin & Tishue [97] performed a detailed comparison of the evolution of preheating between cosmological perturbation theory, rigid FLRW and NR. They used the BSSN formalism with the modified version of the 1+log gauge condition in Eqn. (40). To translate between the results of the three schemes they prepared a useful appendix expressing their Bardeen potential diagnostics in terms of BSSN quantities. They used the variance of fields as a measure of the inhomogeneities, e.g. tracking the variance of the lapse to look for signs of BH formation. Although they found a disagreement with perturbation theory, they were not able to conclude that this led to the formation of BHs.  \nWhen the decay of the inflaton field happens very rapidly, the process is highly non-perturbative and can feature resonances that lead to the fragmentation of the homogeneous condensate into lumps, which then evolve to form oscillons: massive, localized, metastable configurations of scalar field [206], see Fig. 15. Many inflationary potentials can be self-resonant, where resonance generates quanta of the inflaton field itself more efficiently than particles coupled to the inflaton. In a class of well-motivated self-resonant models, the universe may become dominated by oscillons, delaying thermalization. Several questions about oscillons have been explored using NR simulations: (1) Do non-relativistic (and/or nonbackreacting) simulation results still apply when these are taken into account; (2) what is the spectrum of GWs emitted by these objects? (3) can they can collapse to black holes?  \nKou et al. [207], performed NR simulations of oscillon preheating starting from a standard spectrum of the initial vacuum fluctuations ϕ and ˙ ϕ = 0, with an ultraviolet cutoff to ensure convergence. They compared their results using the BSSN formulation with the adapted 1+log gauge conditions in Eqn. (40) implemented within CosmoGRaPH [208] to traditional schemes that use a rigid FLRW background. To quantify whether oscillons dominate the universe, the authors computed the energy fraction contained in oscillons as the fraction of energy in regions where the energy density is greater than twice the average  \nf = ∫ ρ> 2 ⟨ ρ ⟩ dV ρ ∫ dV ρ (59)  \nThey found that for some parameters of the models studied, significant deviations from FLRW were observed (see the left panel of Fig. 16). These were interpreted as GR corrections producing more numerous and denser oscillons, although they could also be a consequence of the gauge and slicing  \n<!-- image -->  \nFig. 16 : (Left:) Comparison of the time evolution of the fractions of energy contained in oscillons between the full GR and FLRW schemes, from Kou et al. [207]. (Right:) Compactness C = GM/R of oscillons formed during inflationary preheating for different scales of the potential µ and initial amplitude of the fluctuations ⟨ δϕ 2 ⟩ , from Aurrekoetxea et al. [73].  \n<!-- image -->  \ndependency of the quantity f . The authors also suggested that these oscillons could collapse into BHs. However, this conclusion was drawn from studying the evolution of a single box-size mode with periodic boundary conditions, a setup that could favour the formation of BHs as it does not permit the oscillons to radiate.  \n̸  \nAurrekoetxea et al. [73] studied preheating with the simplest and most Planck-consistent initial conditions in both the field and velocity at the end of inflation ( ˙ ϕ = 0). It was found that the assumption ˙ ϕ = 0 (mainly the lack of a homogeneous mode, not inhomogeneities) can impact the strength of the resonance and the growth of overdensities. They focused on studying the compactness of the oscillons  \nC = GM R . (60)  \nTheir main result is shown in Fig. 16 (right panel), where the compactness of the oscillons is plotted against µ for different initial amplitude of the fluctuations ⟨ δϕ 2 ⟩ . Here µ is the characteristic scale of the inflationary potential and parameterises whether the model is small field ( µ ≪ M Pl ) or large field ( µ ≈ M Pl ). Their simulations showed that the compactness in bounded by a maximum, which is not sufficient for the oscillons to collapse into BHs. The authors also observed differences in the evolution of overdensities when comparing their results to those obtained using a rigid FLRW background scheme.  \nKou et al. [98] extended their simulations to extract the gravitational waves emitted during and after the formation of these oscillons. In their work, they compared the extracted density ρ GW using a rigid FLRW scheme to NR simulations in different gauges and concluded that the extracted GWs are highly sensitive to the choice of gauge. In particular, the synchronous gauge was found to be a very unreliable gauge condition to be used to extract GW spectra, as it contained redundant gauge modes that artificially enhanced the GW power spectrum. The authors argued that both 1+log and the radiation gauges can provide accurate GW predictions if their parameters are suitably chosen, but further work would be needed to connect these gauge dependent results to observational quantities.  \nAnother example of a preheating process is gauge preheating, when the energy of the inflaton is transferred into gauge fields. This can be sourced via a coupling α g ϕF µν ˜ F µν in the Lagrangian, where α g parameterizes the strength of the coupling and the range of momenta for which the process features a tachyonic instability. This instability allows the full vacuum energy to be transferred within a single oscillation of the inflaton field, and the strongest effects occur when the inflaton rolls the fastest ( ˙ ϕ is the largest) - typically at the end of inflation. In Adshead et al. [102], the authors showed that linearized gravity is quickly violated during this process. However, despite reaching large density contrasts δρ/ρ ≈ 30, PBHs did not form, indicating that the internal pressure of the gauge fields was playing an important role. Using a similar approach to Aurrekoetxea et al. [73] but based on the over-mass δM enclosed in a region of radius R , they computed the compactness of the overdensities formed and found a maximum value of C ≈ 0 . 1. The authors compared the extracted GW spectrum from FLRW and BSSN  \nFig. 17 : Present-day gravitational-wave spectra extracted from simulations with different couplings α g . The figure, from Adshead et al. [102], compares the results obtained from FLRW and full GR simulations.  \n<!-- image -->  \nlattice simulations finding remarkable agreement and suggesting that non-linear gravity at most had a O (1) effect, see Fig. 17.', '4.2 Phase transitions and topological defects': 'As the universe expands and cools down, it can undergo phase transitions that lead to the formation of exotic objects. Among the most studied relics are topological defects [209-211] - structures formed when the universe experienced symmetry breaking. The specific details of these structures depend on the topology of the vacuum. For example, if the vacuum breaks a parity symmetry, two-dimensional defects known as domain walls can form. We previously encountered such relics in Sec. 3.4, where we discussed the observational signatures in the CMB from inflationary bubble collisions. NR has also been used to study the evolution of individual bubbles during vacuum decay, such as the 1+1D work by Blanco-Pillado et al. [212], and Giombi et al. [213]. Recently, Lin et al. [214, 215] have performed 3+1D NR simulations to investigate how the evolution of inhomogeneities in the early universe can populate the landscape of multiple vacua. The dynamics of the formed bubbles may serve as seeds for supermassive black holes, as studied by Huang et al. [216], and enhance the formation of black hole binaries, as suggested by Yuwen et al. [217].  \nIf the transition breaks a U(1) symmetry, line-like defects known as cosmic strings are formed [218, 219]. The vast separation of scales involved in these phenomena already poses a challenge to their study in a rigid expanding background, as the objects quickly fall below the grid resolution. In absence of AMR techniques, it is common to use the fat-string approximation, which introduces an artificial scale that increases the width of the string over time. While this approximation can capture the large-scale dynamics of cosmic strings, it may not accurately represent the small-scale behaviour. There is a need for more sophisticated simulation techniques to resolve these objects across different scales.  \nThe dynamics of topological defects sources source gravitational waves. If these are weak, we might observe the incoherent sum of all the events that are too weak to be directly detected - the stochastic GW background. However, local dynamics might also source strong gravitational-wave signatures that can be directly detected. In both cases, an NR treatment may be necessary: either to validate the predictions of a weak-field treatment of the source, or to predict the detailed signatures to be expected. However, much of the current research focuses on the linearised treatment of these gravitational-wave signatures, and to the best of our knowledge, there is no simulation that fully incorporates NR in a cosmological context -  \nhence the brevity of this section. However, asymptotically flat NR simulations have been performed to understand the GR dynamics of these objects and construct gravitational-wave templates [220-224].', '4.3 Primordial black holes': 'Gravitational collapse of matter has long been a target of NR simulations, in particular the phenomenon of critical collapse uncovered by Choptuik [225] (see [226] for a review). The formation of primordial black holes (PBHs) in the early universe is a case of gravitational collapse that is of particular cosmological interest. PBHs can form through various processes in the early universe, but the primary mechanism considered is the collapse of large primordial fluctuations generated during or at the end of inflation. This is in contrast to ordinary BHs that form as the endpoint of stellar evolution. If sufficiently large fluctuations re-enter the cosmological horizon during the radiation era, they can collapse to form BHs with a wide range of masses. PBHs could be a candidate for dark matter or provide seeds for supermassive BHs. Important questions for NR include quantifying the size and scale of density perturbations that lead to the formation of a BH, and characterising their subsequent evolution, e.g. whether they accrete further matter.  \nAn overdense region is parameterised with respect to the background density by δρ ( R ) ≡ ρ/ ¯ ρ -1. The amplitude of overdensities is commonly defined by integrating this quantity  \nδ ≡ 1 V ∫ R 0 0 4 πδρR 2 dR (61)  \nwhere V = 4 πR 3 0 / 3 and R 0 is the radius of the overdensity. Some early analytical estimates by Carr & Hawking [227, 228] suggested a threshold amplitude for collapse as δ c ≈ 1 / 3, which is significantly above the amplitude of the scale-invariant density fluctuations generated during inflation that we measure on larger (CMB) scales and so requires an enhancement mechanism. Even in such scenarios, PBHs are produced from the high-amplitude tail of the probability distribution of the curvature perturbations. The threshold amplitude δ c therefore needs to be accurately determined, since even a small error in its value can lead to a significant change in the predicted number of PBHs. NR simulations have been used to refine this estimate, and set bounds on the possible production of PBHs in different models.  \nOne of the most common approaches to simulate the collapse of a spherically symmetric overdensity is using the approaches developed by Misner & Sharp [229] and May & White [230]. The spherically symmetric line element is parameterised as  \nds 2 = A 2 ( t, r ) dt 2 + B 2 ( t, r ) dr 2 + R ( t, r ) 2 d Ω 2 , (62)  \nwhere d Ω 2 = dθ 2 +sin 2 θdϕ 2 . The operators D t ≡ A -1 ∂ t and D r ≡ B -1 ∂ r , which represent derivatives with respect to proper time and radial proper distance, are introduced and applied to the areal radius function R ( t, r ) to define the quantities  \nU ≡ D t R = 1 A ∂ t R, Γ ≡ D r R = 1 B ∂ r R. (63)  \nThe quantity U measures the velocity of the fluid with respect to the centre of the coordinates and Γ gives a measure of the spatial curvature. These are related to the Misner-Sharp mass M as  \nΓ 2 = 1 + U 2 -2 M r , (64)  \nwhere M = 4 π ∫ ρR 2 dR . The system of equations for relativistic hydrodynamics is closed after specifying an equation of state that relates the pressure and energy density of the fluid. However, one known problem of the Misner-Sharp formalism is the appearance of coordinate singularities soon after the formation of BHs, which prevents any further evolution.  \nThe first numerical study of PBH formation was conducted by Nadezhin et al. [231]. Using a similar code to that of May & White [230] and incorporating excision to avoid the singularity, the authors explored how the formation process depended on the amplitude and profile of a deviation from FLRW. They established a local criterion for BH formation, determining the proportion of matter collapsing into the PBH and defining its boundary in terms of the visibility horizon . Their results indicated that pressure is a significant impediment to PBH formation, so that when a PBH forms, it is much smaller than the  \n<!-- image -->  \nFig. 18 : (Left:) Time evolution of an overcritical Gaussian perturbation that forms a PBH, simulated using the the Hernandez-Misner formulation. Figure from Niemeyer & Jedamzik [51]. (Right:) Figure from Musco et al. [236] illustrating the evolution of the mass-energy during PBH formation using an explicit Lagrangian hydrodynamics code, in which they found δ c ≈ 0 . 4.  \n<!-- image -->  \ncosmological horizon, and the accretion afterwards is negligible. Novikov & Polnarev [232] investigated the impact of varying the equation of state of the fluid. Softer equations of state resulted in similar behaviours but with more massive PBHs that formed earlier, whereas for hard equations of state with P ≈ ρ , they observed a discontinuity in the flow near the PBH that halted the accretion of matter. Pressure gradients caused the outer layer of matter to be ejected to infinity, and so they concluded that catastrophic accretion (where the PBH grows at the same rate as the cosmological horizon) was impossible. Bicknell & Henriksen [233] used a method of hydrodynamics characteristics to avoid the BH singularity. They showed that regions of overdensity considerably smaller than δρ/ρ ≈ 1 could collapse to PBHs with mass greater than or of the order of the horizon mass.  \nTo avoid the singularity, the Hernandez-Misner [234] version of the Misner-Sharp equations is often used, in which the time coordinate t is replaced by the outgoing null coordinate u . A useful comparison between the Misner-Sharp and the Hernandez-Misner system of equations can be found in Table 1 of Baumgarte et al. [235] or section 2 of Musco et al. [236]. A drawback of this approach is that specifying initial data on a null hypersurface is not trivial, but methods to do this have been developed [235, 237]. Niemeyer & Jedamzik [51] used this approach to investigate the collapse of horizon-sized density fluctuations to PBHs (see Fig. 18). They found that the growth of the PBH mass due to accretion was insignificant, and identified the threshold for BH formation to be δ c ≈ 0 . 7, higher than the commonly used value δ c ≈ 1 / 3. Shibata & Sasaki [52] studied the problem using constant mean curvature timeslicing conditions in 1+1D, and suggested that these differences were due to the presence of an unphysical decaying mode in the initial data. Further work by Green et al. [238], Hawke & Stewart [239] and Musco, Miller & Rezzola [236, 237, 240] confirmed this, finding that the effective δ should be in the range 0 . 3 < δ c < 0 . 5 (see Fig. 18).  \nThese works also considered whether the critical scaling found by Choptuik [225] in asymptotically flat spacetimes and massless scalars  \nM bh = K ( δ -δ c ) γ (65)  \nholds for fluids and in an asymptotically FLRW spacetime, which is the relevant case for PBH formation [241]. Niemeyer & Jedamzik [51] demonstrated that near-critical behaviour applies to PBH formation, extending the results of Evans & Coleman [242] to an asymptotic FLRW background. Their simulations yielded the critical exponent γ ≈ 0 . 36 and K ranging from 2 . 7 to 12 depending on the shape of the initial overdensity. Hawke & Stewart found consistent results, but noted that the scaling relation did not hold for sufficiently small ( δ -δ c ), reaching a minimum PBH mass M bh ≈ 10 -4 M h . Musco et al. [236] also confirmed the scaling law relation and found that a positive cosmological constant leads to lower values of γ and higher critical δ c .  \nWork has continued using the Misner-Sharp and Hernandez-Misner equations to study PBH formation in 1+1D in a variety of contexts [243-259]. The 3+1D NR techniques that are now standard in the simulation of BBH mergers have not yet been extensively used by the PBH community, with the exception of the early work by Shibata & Sasaki [260], and simulations performed with COSMOS [261-263] (based on the NR code SACRA [264]) and GRChombo [265] NR codes. These techniques have been used to go beyond 1+1D simulations and study the effect that non-spherical perturbations may have on the results.  \nUsing spheroidal super-horizon perturbations, Yoo et al. [266] studied the effect of ellipticity on the threshold of PBH formation, finding that it is negligible in standard situations during the radiationdominated era. Yoo later extended the work [267] by adding a tidal torque to the non-spherical profile of the initial curvature perturbation as a mechanism to generate angular momentum and form a spinning remnant, but the final spin remained consistent with zero, in agreement with expectations for the radiation-dominated era. However, in the case of a softer equation of state, small deviations from spherical symmetry may evolve differently due to the absence of pressure gradient forces. De Jong et al. [74, 268] used a simple model of two scalar fields to study the mass and spin of PBHs from the collapse of perturbations with angular momentum during an early matter-dominated era P ≈ 0. They found that although the initial mass and spin at formation can be of order O (10%), the accretion of background matter post-collapse dilutes the dimensionless spin of the PBH. Hence, unless the matter era is short, they concluded that the final spin of the PBHs formed will be negligible.  \nThere appears to be significant scope for these studies to be extended to more general cases, using relativistic hydrodynamics and modern 3+1D NR methods.', '5 Revelation: The late universe': "Observations of large-scale structure in the cosmology of the late universe - by which we broadly mean cosmological evolution past the epoch of matter-radiation equality - provide a key probe of cosmological parameters and fundamental physics. Observations of the CMB [1, 269], galaxy clustering and matter distribution [16] and galaxy lensing [270] constrain the parameters of the ΛCDM model of cosmology, which is presently considered the standard description of our universe.  \nMuch of the interpretative power of these observations come from their statistics , e.g. their 2-point and higher correlation functions. While some predictions can be made analytically, the gold standard results are derived from numerical simulations of large-scale structure using N-body/fluid codes (for a review see [271]). These simulations encompass scales on which the large-scale (Hubble-Lema ˆ i tre) space-time curvature plays a role, but usually make several simplifying assumptions compared to a full NR evolution. First, it is assumed that the matter component is non-relativistic and (for the dark matter) collisionless. In addition, since on the largest scales the universe is approximately homogeneous it is assumed that the gravitational fields of perturbations are weak. This means that, locally at least, gravitational forces are captured by Newtonian gravity, specifically the Poisson equation,  \n∇ 2 Φ = 4 πGρ , (66)  \nwhere Φ is the local gravitational field sourced by the matter density ρ . Dark matter is represented by massive point particles that are evolved using Newton's 2nd Law of motion a = ∇ Φ where a is the acceleration. This so-called Poisson-Vlasov limit , is employed by most of the main N-body simulation codes e.g. GADGET-4 [272] and RAMSES [4]. Since the effects of expansion are important at scales comparable to the Hubble scale H -1 even in the Poisson-Vlasov limit, most N-body codes capture this effect by allowing the background space to homogeneously expand, following the Friedman equations of some pre-specified ΛCDM cosmology.  \nThere are two related motivations for going beyond these assumptions in late-time cosmology with NR. First, to quantify the effects of the backreaction of the (possibly O (1)) small-scale inhomogeneities on the large-scale evolution. This was initially related to the idea that small-scale inhomogeneities might somehow manifest on the largest scales as a smooth negative pressure fluid (hence mimicking a cosmological constant) [273, 274]. Part of this question involves understanding how such inhomogeneities will impact on actual observations by comoving (or in the more general case geodesic) observers. Second, higher order correlation functions beyond the power spectrum can be sourced by non-linear evolution, thus in principle any search for non-Gaussianities [275-278] in the large-scale structure data must account for non-linear effects. Even if non-linear effects are too small to be detected, the increasing size and precision of large-scale structure datasets, especially in the upcoming Stage 4 class surveys, means that all  \n<!-- image -->  \n<!-- image -->  \nFig. 19 : (Left:) Profile of the inhomogeneous matter density ratio ρ/ ¯ ρ on the y = z plane at the start of the simulation and at a later time, from Bentivegna & Bruni [84]. (Right:) The ratio of the proper length to the initial length of six sets of paths in the simulation, compared with the FLRW estimate, from shorter paths (purple) to longer (red) in spectral order. Shorter paths deviate more from the FLRW expectation, up to a maximum of ∼ 5% over one e-fold. From Giblin et al. [79].  \n<!-- image -->  \npossible systematics must be accounted for [279]. This makes an investigation of large-scale structure evolution with NR highly timely.  \nThere exist several codes that are half-way houses between the Poisson-Vlasov approach and full NR, which are worth highlighting since they capture many of the potential relativistic effects with fewer overheads than a full NR simulation. For example, in gevolution [280], the Poisson equation (which arises from the t -t component of the Einstein equations in the non-relativistic weak-field limit), is extended by the full t -t and i -j Einstein equations up to quadratic order, to better capture the dynamics of the gravitational degrees of freedom. The first equation is an elliptic equation that captures the dynamics of the Newtonian potential, while the second captures the tensor degrees of freedom as a hyperbolic equation. Another code is GRAMSES [281, 282], in which the metric degrees of freedom are obtained by solving the momentum and Hamiltonian constraints in-between time steps (instead of using the Poisson equation). This stratagem captures the gravitational degrees of freedom in the non-relativistic limit. In both codes, the matter sector is modelled by collisionless particles that move according to their corresponding geodesic equations.", '5.1 Backreaction and averaging in non FLRW spacetimes': "The first systematic investigations using NR to study late-time cosmological evolution were done using the CosmoGRaPH code of Giblin, Mertens & Starkman [79, 208], and in parallel by Bentivegna & Bruni [84] using the McLachlan evolution thorn of Einstein Toolkit suite of tools [283, 284]. See Fig. 19.  \nIn these works, matter was modelled as a pressureless fluid with density parameter ρ , and evolved via the continuity equation - that is, matter was assumed to flow along geodesics, with zero initial coordinate velocity, such that the momentum flux S i = 0 identically. Both works successfully used the synchronous gauge ( α = 1 , β i = 0), which does not completely remove issues of gauge dependence in inhomogeneous spacetimes but means that the coordinate observers track geodesics. Therefore in the limit of small perturbations one should recover the FLRW slicing in the comoving gauge. In a later paper Giblin et al. [285] used ray-tracing of null geodesics in the same simulations to properly quantify the effects of inhomogeneities on observables and to construct an average Hubble diagram , and found consistent results.  \nIn order to match numerical simulations to the FLRW quantities, the authors defined volume averaged quantities ⟨ X ⟩ S over the spatial hyperslices S  \n⟨ X ⟩ S ≡ ∫ S X √ γd 3 x ∫ S √ γd 3 x , (67)  \nwhere the normalisation factor in the denominator is simply the spatial volume, from which we can obtain the approximate scale factor-like average  \na 3 S ≡ ∫ S √ γd 3 x . (68)  \nIn Giblin et al. [79] it was shown that the averaged quantities closely tracked those of an equivalent FLRW spacetime. When perturbations (in the form of a peaked power spectrum) with variance ⟨ δρ ⟩ / ⟨ ρ ⟩ ≈ 0 . 3  \nwere introduced, they found that the deviation from FLRW, captured in the variance of the local expansion parameter ⟨ δK ⟩ / ⟨ K ⟩ was of order 0 . 3% suggesting that even in the presence of large perturbations, spacetime tracks FLRW to a good approximation.  \nIn Bentivegna & Bruni [84, 286], the authors quantified the deviation from FLRW by computing the correction terms to the Friedman equation. Taking the 2nd time derivative of Eqn. (68), yields an equation of motion that is analogous to the 2nd Friedman equation  \na S a S = -4 π 3 ⟨ ρ ⟩ a -3 S + Q S 3 , (69)  \nwhere the correction term  \nS = 2 3 ( ⟨ K 2 ⟩ - ⟨ K ⟩ 2 -2 ⟨ A 2 ⟩ ) , (70)  \nmeasures the back-reaction from non-linear evolution. They found that for small values of initial perturbations, Q D ∝ a -1 D , which agrees with 2nd order perturbation theory calculations [287]. They checked the validity of the standard semi-analytical model of collapse, the spherical top-hat model, in which a top-hat shaped overdensity is allowed to evolve under global expansion in a matter-dominated universe. The competition between infall and expansion leads to a moment of turn around when the gravitational collapse overcomes the local expansion. The critical density at turnaround δ OD ≡ δρ/ρ at the turnaround time t TA . They found that δ OD differed to the semi-analytic estimate of δ OD ≈ 1 . 06, but a later work by Munoz & Bruni [288, 289] resolved this difference, finding δ OD = 1 . 055 in good agreement with the semi-analytic results, thus validating the standard collapse model.  \nThe main results of these two works were subsequently confirmed by Macpherson et al. [290], who also used the Einstein Toolkit but incorporated the GRHydro thorn for more general matter and developed a new FLRWSolver thorn to generate the initial data and diagnostics. In this and later works [291, 292] Macpherson et al. simulated late-time structure formation in a matter-dominated universe, using a realistic matter power spectrum of perturbations derived from observations of the CMB as initial conditions. Due to limited resolution, they chose to simulate cosmological scales ranging from k ∼ 10 -2 Mpc -1 to k ∼ 1 Mpc -1 . By implementing an averaging scheme which depends on the scale r D , they found that while large-scale spatial homogeneity and isotropy is maintained, at scales r D ≤ 180 h -1 Mpc, which correspond to local galaxy survey scales, deviation from homogeneity could be up to 6%. Nonetheless, this deviation only affects the measured value of the Hubble constant at the sub-percent level, and so is not sufficient on its own to resolve the current tension in its value, although in Macpherson [293] the value is found to be more significant for the case of observers located close to inhomogeneities. A study that added tensor perturbations was carried out by Wang [294].  \nThe codes used in these works undertook a comparison exercise [295] together with gevolution and GRAMSES . All codes ran a simulation with identical initial conditions and similar resolutions and compared the paths of two null, spatially separated trajectories, launched from the initial hyperslice. As the null trajectory is sensitive the entire spacetime metric, this serves as a good test of equivalence in relativistic cases. All codes showed excellent agreement (within 0 . 03%).  \nA direct comparison between NR and N-body simulations was carried out by East et al. [296]. In this work, a full hydrodynamical treatment was used to model the collisionless dark matter [297], while the metric variables were evolved using the GHC formulation [90]. Meanwhile, GADGET-2 [3] was used to generate results for a non-relativistic N-body simulation for comparison to the NR results. Instead of computing deviations from linear theory using averaged quantities, the authors chose a set of fiducial observers that were co-moving with the matter, and tracked the evolution of δρ over their worldlines. They compared their results to both the Poisson-Vlasov results and expectations from linear perturbation theory. Furthermore, they solved for the geodesics of null rays, and therefore allowing them to compute the cosmological distance as a function of redshift for free-falling observers located along the photon paths. They found that, as long as the Newtonian potential remains in the weak-field limit, the difference between the NR simulations to both the N-body simulations and linear theory expectation is at the sub-percent level, see Fig. 20.  \nWhere conservative hydrodynamical equations are used to model the dark matter, as in East et al. [296], direct comparison between the fluid picture employed in the NR simulations and the particle picture in the N-body simulations become unreliable once shell-crossing occurs, as the virialisation process cannot be captured in the fluid description. To overcome this, one can introduce the Einstein-Vlasov formalism, in which matter is modelled as a distribution of N collisionless point-like particles. In this formalism,  \n<!-- image -->  \nFig. 20 : (Left:) The density contrast as measured by an observer comoving with matter at the point of maximum overdensity (OD) or underdensity (UD) for two cases and different initial velocity profiles. The GR (solid lines) and Newtonian (dotted lines) calculations closely track each other. Using fluid simulations in GR, from East et al. [296]. (Right:) The same results but for all different treatments including the Einstein-Vlassov system (blue and green lines), again showing good agreement between the Newtonian and GR cases. From East et al. [298].  \n<!-- image -->  \nthe stress-tensor is given by the sum of all the energy of all particles labelled by i  \nT µν = i = N ∑ i =1 ϵ i l µ i l ν i , (71)  \nwhere ϵ i is the energy density of the particle, and l µ i is the tangent to geodesic curve of the particle x µ i ( τ ). For particles of mass m , l µ i l iµ = -m 2 . The particles follow the geodesic equation as usual,  \nl µ i ∇ µ l ν i = 0 . (72)  \nAn averaging prescription is then implemented to coarse-grain Eqn. (71) into an average value per grid point. This was first implemented by Pretorius & East [299] to study the collapse of collisionless massless particles and extended to massive particles in a later work [298]. In parallel, Giblin et al. [300] implemented this formalism into the CosmoGRaPH framework for collisionless massive particles, and more recently by Magnall et al. [301] using the phantom SPH code.  \nThe Einstein-Vlasov framework allows one to directly compare along the trajectories of particles beyond shell crossing. As in the earlier work with fluids [296], these works found that the Einstein-Vlasov formalism agrees with the Poisson-Vlasov formalism up to sub-percent level, as long as the Newtonian potential does not exceed the weak-field approximation | Ψ N | ≤ 0 . 1. In East et al. [298], setting up identical initial conditions in both formalisms allowed them to explicitly compare the trajectories, and they found that as long as the evolution remained non-relativistic (in the Einstein-Vlasov picture), they were in excellent agreement with each other. Interestingly, they found that in the non-linear limit, the collapse of overdensities occurred faster in the Newtonian picture when compared to the GR one, although even then differences remained small, see Fig. 20.  \nThe use of NR to probe the dynamics of expanding universes and the impact of backreaction in inhomogeneous cases has also been tackled from the perspective of a lattice of BHs. Roughly speaking, if the scale of the BHs is much shorter than the curvature scale of the universe M BH ≪ H -1 , then a lattice of static BHs mimics (at the largest scales) that of cold dark matter (see e.g. [302]), which can drive cosmic expansion. Thus such a spacetime is an interesting toy model , with which one can study the physics of dark matter and the role of averaging in FLRW, at least on the largest scales. This was first proposed by Linquist and Wheeler [303], and then numerically simulated in several works [304310], see Bentivegna et al. [311] for a review. These simulations have studied effects such as the dressing (differences from the naive superposition of constituents) of mass and other physical parameters, and the biases in optical measurements that inhomogeneous gravitational fields can produce.  \nOther notable work on averaging focusing on a confirmation of Isaacson's formula (that gravitational waves behave as a radiation fluid on average) was performed by Ikeda et al. [312].", '5.2 Cosmological observables in the non-linear regime': "A general theme of the results described in the previous works in this section is that the PoissonVlasov formalism is remarkably accurate in capturing most of the relevant gravitational effects in latetime cosmological evolution, even in the non-linear regime . Significant deviations from full NR do not occur until deep in the strong gravity regime - expected at very short wavelength scales when local regions decouple fully from the Hubble flow. Indeed, these results indicate that scale-separation is robust in cosmological evolution - non-linearity of GR does not induce UV-IR couplings where UV modes significant affect large-scale cosmological evolution.  \nNevertheless, as mentioned earlier, NR still has a role to play both in terms of understanding the systematics of N -body predictions and as a useful cross-check. The first such cross-check was undertaken by Giblin et al. [313] where they computed the weak lensing power spectrum using NR. Weak lensing power computations often rely on approximations [314]. In particular, line-of-sight integrals often rely on the so-called Born approximation , where the light ray is assumed to propagate along the background geodesic and is not perturbed by higher order gravitational effects. Attempts to incorporate post-Born effects has led to disagreement 13 [317, 318]. Giblin et al. [313] attempted to resolve this issue by computing the weak lensing power spectrum by ray-tracing a universe evolved using NR. They found that the approximation may lead to percent level error at angular scales of l = 20 ∼ 30 when comparing full NR to perturbation theory computations, potentially signalling important systematics not captured in higher order perturbation theory. More recent efforts to develop tools to probe cosmological distances using ray tracing have been carried out [319, 320], and support the conclusion that effects may be at percent level for scales above 100 Mpc.  \nMacpherson & Heinesen [321] evolved an anisotropic cosmological spacetime, and then computed observables such as the Hubble parameter and its higher moments with different levels of sky coverage. They showed that there was up to percent level cosmic variance in the measurements. For example, there is a 7% probability that an observer can measure up to a 3% difference in the Hubble parameter of a typical Einstein-de-Sitter universe, see Fig. 21. In related works, Macpherson [293] and Williams et al. [322] have studied the impact of nearby inhomogeneities and voids on observers' inferences of the Hubble constant, where they found that the voids are in general negatively curved and void statistics such as size distribution roughly follows that of N -body simulations. Koksbang et al. [323] extended this by computing the red-shift drift - where the redshift of an object evolves through time due to cosmic expansion. They showed that NR predicts that at redshift z = 0 . 1, fluctuations in the redshift drift due to inhomogeneities are at the ∼ 10% level, and so potentially observable in the next generation spectroscopic sky surveys.  \nAnother use of NR has been to probe the question of cosmological spatial curvature Ω k - how spatially flat is the universe - in an inhomogeneous universe, as in Tian et al. [324]. Since curvature perturbations violate FLRW symmetry, averaged quantities over a spatial hyperslice do not necessarily correspond to the true local spatial curvature. The authors found that these inhomogeneities introduce an effective Ω k of the order of the perturbations, in general agreement with analytic arguments. This has recently been extended by Oestreicher & Koksbang [325] who found more significant curvature perturbations of order 10%.", '6 Closing remarks': 'Cosmologists seem to have struck lucky. They profit from the simplicity of the FLRW model, due to the approximate homogeneity and isotropy of our observable universe. Beyond this realm, where the light of the CMB peters out and gravity may be the only messenger, there is no guarantee that the FLRW model still holds. To explore this regime we are forced to rely on new methods, and NR provides a particularly useful tool where dynamical strong gravity effects are at play. The use of NR introduces conceptual complications, challenges to implementation, and a potentially high cost per simulation. But used efficiently it can (and already has) shed light on parts of the Universe that are dark to us.  \nIn this review, our aim is twofold. First, we provide a starting point for both cosmologists and numerical relativists who are interested in using NR in the context of cosmology. While much of the specialised NR knowledge developed for BH spacetimes carries over to cosmology, the mapping is not perfect and there is an opportunity for new techniques to be developed. Finding physically well-motivated  \nFig. 21 : Sky-maps of the effective Hubble, deceleration, curvature, and jerk parameters (top-left to bottom-right, respectively) for an observer. Each parameter (with the exception of the curvature) is normalised by its respective EdS value. From Macpherson & Heinesen [321].  \n<!-- image -->  \nconstraint-obeying initial conditions and developing gauge conditions adapted to cosmological spacetimes far from the FLRW universe are two good examples of areas that deserve further study.  \nSecond, we give a comprehensive overview of cosmological simulations that have used NR, covering a range of work: from tentative steps to explore the unknown darkness beyond the observable universe, to important work checking GR backreaction effects in numerical simulations of the late universe. There is much more work to be done, and many areas where we have barely scratched the surface - exploring cosmology with modified gravity, understanding string cosmology, extending gravitating cosmic strings simulations to networks, to name a few.  \nWe hope that this Living Review will encourage cosmologists and numerical relativists alike to break new ground, and we look forward to expanding this Review in the future with the fruits of their labour.', '7 Summary table and open questions': 'A summary table of the work that has been completed in the various areas covered by the review is provided in Tab. 1 for reference. We also provide here an (incomplete) list of open questions and outstanding directions that have been highlighted by the review:', 'NR techniques for cosmology': "- 1. What are the optimum boundary and gauge conditions for cosmological simulations, and how do they depend on the physics being studied?\n- 2. How can we access a wider range of physical scales or make approximations to allow us to study GR effects across scales, e.g. in string networks?\n- 3. What are the best diagnostics for inhomogeneous spacetimes? Can spatial averages be trusted or do we always need to use geodesic observers? If so, how should they be positioned to sample the spacetime?\n- 4. How can we accurately quantify GWs or other stochastic background observables in cosmological spacetimes?  \nTable 1 : Summary of the work in this review.  \n- 5. How can we construct initial data that is physically 'correct', as well as constraint satsfying? To what extent do arbitrary choices like conformal flatness limit the physics we can study? Can we go beyond these?\n- 6. Can we develop open source tools that are general purpose for different cosmological spacetimes, e.g. a general initial condition solver for 3+1D periodic spacetimes?", 'Genesis: the early universe': '- 1. Does the BKL conjecture continue to hold for asymptotically flat spacetimes?\n- 2. What is the role of periodic boundary conditions in findings about the robustness of inflation simulations? Do they provide a bias towards inflation?\n- 3. What is the impact of considering multi-field models on the robustness of inflation?\n- 4. How can the smooth contraction phase connect to the bounce in Ekpyrosis scenarios?\n- 5. What are the observable consequences of either inflation or bouncing scenarios in the CMB or other observations like GW backgrounds?\n- 6. Do modifications to gravity at higher energies change these behaviours? Can we connect to quantum gravity/string theory predictions?', 'Exodus: the transitional universe': '- 1. How does gravity and AMR impact the evolution of a network of early universe relics? e.g. oscillons, topological defects.\n- 2. Are there any (p)reheating scenarios that generate significant amounts of PBHs with small amplitude fluctuations?  \n- 3. How can we best quantify the stochastic GW background produced from early universe phenomena like reheating?\n- 4. What are the effects of non-spherical perturbations (i.e. going to 3+1D simulations) in the formation of PBHs?\n- 5. Can we make quantitative predictions that connect observational data to parameters of the fundamental theory?', 'References': "- [1] Aghanim, N., et al. : Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. 641 , 6 (2020) https://doi.org/10.1051/0004-6361/201833910 arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)]\n- [2] Baumann, D.: Cosmology. Cambridge University Press, Cambridge (2022). https://doi.org/10. 1017/9781108937092\n- [3] Springel, V.: The Cosmological simulation code GADGET-2. Mon. Not. Roy. Astron. Soc. 364 , 1105-1134 (2005) https://doi.org/10.1111/j.1365-2966.2005.09655.x arXiv:astro-ph/0505010\n- [4] Teyssier, R.: Cosmological hydrodynamics with adaptive mesh refinement: a new high resolution code called RAMSES. Astron. Astrophys. 385 , 337-364 (2002) https://doi.org/10.1051/0004-6361: 20011817 arXiv:astro-ph/0111367\n- [5] Anninos, P.: Computational cosmology: From the early universe to the large scale structure. Living Rev. Rel. 4 , 2 (2001) https://doi.org/10.12942/lrr-2001-2 arXiv:gr-qc/0108083\n- [6] Berger, B.K.: Numerical approaches to space-time singularities. Living Rev. Rel. 5 , 1 (2002) https: //doi.org/10.12942/lrr-2002-1 arXiv:gr-qc/0201056\n- [7] Andersson, L.: On the relation between mathematical and numerical relativity. Class. Quant. Grav. 23 , 307-318 (2006) https://doi.org/10.1088/0264-9381/23/16/S02 arXiv:gr-qc/0607065\n- [8] Cardoso, V., et al. : NR/HEP: roadmap for the future. Class. Quant. Grav. 29 , 244001 (2012) https://doi.org/10.1088/0264-9381/29/24/244001 arXiv:1201.5118 [hep-th]\n- [9] Lehner, L., Pretorius, F.: Numerical Relativity and Astrophysics. Ann. Rev. Astron. Astrophys. 52 , 661-694 (2014) https://doi.org/10.1146/annurev-astro-081913-040031 arXiv:1405.4840 [astroph.HE]\n- [10] Choptuik, M.W., Lehner, L., Pretorius, F.: Probing Strong Field Gravity Through Numerical Simulations (2015) arXiv:1502.06853 [gr-qc]\n- [11] Garfinkle, D.: Numerical Relativity Beyond Astrophysics. Rept. Prog. Phys. 80 (1), 016901 (2017) https://doi.org/10.1088/0034-4885/80/1/016901 arXiv:1606.02999 [gr-qc]  \n- [12] Coley, A., Lehner, L., Pretorius, F., Wiltshire, D.: Computational Issues in Mathematical Cosmology (2017)\n- [13] Ijjas, A.: Numerical Relativity as a New Tool for Fundamental Cosmology. MDPI Physics 4 (1), 301-314 (2022) https://doi.org/10.3390/physics4010021 arXiv:2201.03752 [gr-qc]\n- [14] Weinberg, S.: Cosmology, (2008)\n- [15] Mukhanov, V.: Physical Foundations of Cosmology. Cambridge University Press, Cambridge (2005). https://doi.org/10.1017/CBO9780511790553\n- [16] Adame, A.G., et al.: DESI 2024 VI: Cosmological Constraints from the Measurements of Baryon Acoustic Oscillations (2024) arXiv:2404.03002 [astro-ph.CO]\n- [17] Riess, A.G., et al. : A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s -1 Mpc -1 Uncertainty from the Hubble Space Telescope and the SH0ES Team. Astrophys. J. Lett. 934 (1), 7 (2022) https://doi.org/10.3847/2041-8213/ac5c5b arXiv:2112.04510 [astro-ph.CO]\n- [18] Freedman, W.L., Madore, B.F., Hoyt, T., Jang, I.S., Beaton, R., Lee, M.G., Monson, A., Neeley, J., Rich, J.: Calibration of the Tip of the Red Giant Branch (TRGB) (2020) https://doi.org/10. 3847/1538-4357/ab7339 arXiv:2002.01550 [astro-ph.GA]\n- [19] Bardeen, J.M.: Gauge Invariant Cosmological Perturbations. Phys. Rev. D 22 , 1882-1905 (1980) https://doi.org/10.1103/PhysRevD.22.1882\n- [20] Ijjas, A., Pretorius, F., Steinhardt, P.J.: Stability and the Gauge Problem in Non-Perturbative Cosmology. JCAP 01 , 015 (2019) https://doi.org/10.1088/1475-7516/2019/01/015 arXiv:1809.07010 [gr-qc]\n- [21] Turner, M.S.: The Road to Precision Cosmology (2022) https://doi.org/10.1146/ annurev-nucl-111119-041046 arXiv:2201.04741 [astro-ph.CO]\n- [22] Maldacena, J.M.: Non-Gaussian features of primordial fluctuations in single field inflationary models. JHEP 05 , 013 (2003) https://doi.org/10.1088/1126-6708/2003/05/013 arXiv:astroph/0210603\n- [23] Weinberg, S.: Quantum contributions to cosmological correlations. Phys. Rev. D 72 , 043514 (2005) https://doi.org/10.1103/PhysRevD.72.043514 arXiv:hep-th/0506236\n- [24] Mukhanov, V.F.: Quantum Theory of Gauge Invariant Cosmological Perturbations. Sov. Phys. JETP 67 , 1297-1302 (1988)\n- [25] Mukhanov, V.F., Feldman, H.A., Brandenberger, R.H.: Theory of cosmological perturbations. Part 1. Classical perturbations. Part 2. Quantum theory of perturbations. Part 3. Extensions. Phys. Rept. 215 , 203-333 (1992) https://doi.org/10.1016/0370-1573(92)90044-Z\n- [26] Gerlach, U.H., Sengupta, U.K.: GAUGE INVARIANT PERTURBATIONS ON MOST GENERAL SPHERICALLY SYMMETRIC SPACE-TIMES. Phys. Rev. D 19 , 2268-2272 (1979) https://doi. org/10.1103/PhysRevD.19.2268\n- [27] Arnowitt, R.L., Deser, S., Misner, C.W.: The Dynamics of general relativity. Gen. Rel. Grav. 40 , 1997-2027 (2008) https://doi.org/10.1007/s10714-008-0661-1 arXiv:gr-qc/0405109\n- [28] Baumgarte, T.W., Shapiro, S.L.: Numerical Relativity: Starting from Scratch. Cambridge University Press, Cambridge (2021). https://doi.org/10.1017/9781108933445\n- [29] Alcubierre, M.: Introduction to 3+1 Numerical Relativity. International Series of Monographs on Physics. Oxford University Press, Oxford (2008). https://doi.org/10.1093/acprof:oso/ 9780199205677.001.0001\n- [30] Baumgarte, T.W., Shapiro, S.L.: Numerical Relativity: Solving Einstein's Equations on  \nthe Computer. Cambridge University Press, Cambridge (2010). https://doi.org/10.1017/ CBO9781139193344  \n- [31] Gourgoulhon, E.: 3+1 formalism and bases of numerical relativity (2007) arXiv:gr-qc/0703035\n- [32] Shibata, M.: Numerical Relativity, (2016). https://doi.org/10.1142/9692\n- [33] Hilditch, D.: An Introduction to Well-posedness and Free-evolution. Int. J. Mod. Phys. A 28 , 1340015 (2013) https://doi.org/10.1142/S0217751X13400150 arXiv:1309.2012 [gr-qc]\n- [34] Shibata, M., Nakamura, T.: Evolution of three-dimensional gravitational waves: Harmonic slicing case. Phys. Rev. D 52 , 5428-5444 (1995) https://doi.org/10.1103/PhysRevD.52.5428\n- [35] Baumgarte, T.W., Shapiro, S.L.: On the numerical integration of Einstein's field equations. Phys. Rev. D 59 , 024007 (1998) https://doi.org/10.1103/PhysRevD.59.024007 arXiv:gr-qc/9810065\n- [36] Nakamura, T., Oohara, K., Kojima, Y.: General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes. Prog. Theor. Phys. Suppl. 90 , 1-218 (1987) https://doi.org/10. 1143/PTPS.90.1\n- [37] Pretorius, F.: Numerical relativity using a generalized harmonic decomposition. Class. Quant. Grav. 22 , 425-452 (2005) https://doi.org/10.1088/0264-9381/22/2/014 arXiv:gr-qc/0407110\n- [38] Pretorius, F.: Evolution of binary black hole spacetimes. Phys. Rev. Lett. 95 , 121101 (2005) https: //doi.org/10.1103/PhysRevLett.95.121101 arXiv:gr-qc/0507014\n- [39] Garfinkle, D.: Harmonic coordinate method for simulating generic singularities. Phys. Rev. D 65 , 044029 (2002) https://doi.org/10.1103/PhysRevD.65.044029 arXiv:gr-qc/0110013\n- [40] Uggla, C., Elst, H., Wainwright, J., Ellis, G.F.R.: The Past attractor in inhomogeneous cosmology. Phys. Rev. D 68 , 103502 (2003) https://doi.org/10.1103/PhysRevD.68.103502 arXiv:grqc/0304002\n- [41] Garfinkle, D.: Numerical simulations of general gravitational singularities. Class. Quant. Grav. 24 , 295-306 (2007) https://doi.org/10.1088/0264-9381/24/12/S19 arXiv:0808.0160 [gr-qc]\n- [42] Estabrook, F.B., Robinson, R.S., Wahlquist, H.D.: Hyperbolic equations for vacuum gravity using special orthonormal frames. Class. Quant. Grav. 14 , 1237-1247 (1997) https://doi.org/10.1088/ 0264-9381/14/5/025 arXiv:gr-qc/9703072\n- [43] Buchman, L.T., Bardeen, J.M.: A Hyperbolic tetrad formulation of the Einstein equations for numerical relativity. Phys. Rev. D 67 , 084017 (2003) https://doi.org/10.1103/PhysRevD.72.049903 arXiv:gr-qc/0301072. [Erratum: Phys.Rev.D 72, 049903 (2005)]\n- [44] Hilditch, D.: Dual Foliation Formulations of General Relativity (2015) arXiv:1509.02071 [gr-qc]\n- [45] Zilhao, M., Cardoso, V., Gualtieri, L., Herdeiro, C., Sperhake, U., Witek, H.: Dynamics of black holes in de Sitter spacetimes. Phys. Rev. D 85 , 104039 (2012) https://doi.org/10.1103/PhysRevD. 85.104039 arXiv:1204.2019 [gr-qc]\n- [46] Daverio, D., Dirian, Y., Mitsou, E.: A numerical relativity scheme for cosmological simulations. Class. Quant. Grav. 34 (23), 237001 (2017) https://doi.org/10.1088/1361-6382/aa9312 arXiv:1611.03437 [gr-qc]\n- [47] Kazdan, J.L., Warner, F.W.: Scalar curvature and conformal deformation of Riemannian structure. Journal of Differential Geometry 10 (1), 113-134 (1975) https://doi.org/10.4310/jdg/1214432678\n- [48] Aurrekoetxea, J.C., Clough, K., Lim, E.A.: CTTK: a new method to solve the initial data constraints in numerical relativity. Class. Quant. Grav. 40 (7), 075003 (2023) https://doi.org/10.1088/ 1361-6382/acb883 arXiv:2207.03125 [gr-qc]  \n- [49] Garfinkle, D., Ijjas, A., Steinhardt, P.J.: Initial conditions problem in cosmological inflation revisited. Phys. Lett. B 843 , 138028 (2023) https://doi.org/10.1016/j.physletb.2023.138028 arXiv:2304.12150 [gr-qc]\n- [50] Coley, A.A., Lim, W.C.: Periodic boundary conditions and G 2 cosmology. Class. Quant. Grav. 41 (1), 015009 (2024) https://doi.org/10.1088/1361-6382/ad0b9f arXiv:2312.03285 [gr-qc]\n- [51] Niemeyer, J.C., Jedamzik, K.: Dynamics of primordial black hole formation. Phys. Rev. D 59 , 124013 (1999) https://doi.org/10.1103/PhysRevD.59.124013 arXiv:astro-ph/9901292\n- [52] Shibata, M., Sasaki, M.: Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity. Phys. Rev. D 60 , 084002 (1999) https://doi.org/10.1103/PhysRevD. 60.084002 arXiv:gr-qc/9905064\n- [53] Corman, M., East, W.E., Ripley, J.L.: Evolution of black holes through a nonsingular cosmological bounce. JCAP 09 , 063 (2022) https://doi.org/10.1088/1475-7516/2022/09/063 arXiv:2206.08466 [gr-qc]\n- [54] Shibata, M., Nakao, K.-i., Nakamura, T., Maeda, K.-i.: Dynamical evolution of gravitational waves in the asymptotically de Sitter space-time. Phys. Rev. D 50 , 708-719 (1994) https://doi.org/10. 1103/PhysRevD.50.708\n- [55] Calabrese, G., Lehner, L., Tiglio, M.: Constraint preserving boundary conditions in numerical relativity. Phys. Rev. D 65 , 104031 (2002) https://doi.org/10.1103/PhysRevD.65.104031 arXiv:grqc/0111003\n- [56] Bona, C., Ledvinka, T., Palenzuela-Luque, C., Zacek, M.: Constraint-preserving boundary conditions in the Z4 numerical relativity formalism. Class. Quant. Grav. 22 , 2615-2634 (2005) https://doi.org/10.1088/0264-9381/22/13/007 arXiv:gr-qc/0411110\n- [57] Sarbach, O., Tiglio, M.: Boundary Conditions for Einstein's Field Equations: Mathematical and Numerical Analysis. J. Hyperbol. Diff. Equat. 2 (04), 839-883 (2005) https://doi.org/10.1142/ S0219891605000634 arXiv:gr-qc/0412115\n- [58] Buchman, L.T., Sarbach, O.C.A.: Towards absorbing outer boundaries in general relativity. Class. Quant. Grav. 23 , 6709-6744 (2006) https://doi.org/10.1088/0264-9381/23/23/007 arXiv:grqc/0608051\n- [59] Rinne, O., Lindblom, L., Scheel, M.A.: Testing outer boundary treatments for the Einstein equations. Class. Quant. Grav. 24 , 4053-4078 (2007) https://doi.org/10.1088/0264-9381/24/16/ 006 arXiv:0704.0782 [gr-qc]\n- [60] Rinne, O., Buchman, L.T., Scheel, M.A., Pfeiffer, H.P.: Implementation of higher-order absorbing boundary conditions for the Einstein equations. Class. Quant. Grav. 26 , 075009 (2009) https: //doi.org/10.1088/0264-9381/26/7/075009 arXiv:0811.3593 [gr-qc]\n- [61] Peterson, C., Gautam, S., Rainho, I., nuales, A., Hilditch, D.: 3D evolution of a semilinear wave model for the Einstein field equations on compactified hyperboloidal slices. Phys. Rev. D 108 (2), 024067 (2023) https://doi.org/10.1103/PhysRevD.108.024067 arXiv:2303.16190 [gr-qc]\n- [62] Gautam, S., nuales, A., Hilditch, D., Bose, S.: Summation by Parts and Truncation Error Matching on Hyperboloidal Slices. Phys. Rev. D 103 (8), 084045 (2021) https://doi.org/10.1103/PhysRevD. 103.084045 arXiv:2101.05038 [gr-qc]\n- [63] nuales, A., Husa, S., Hilditch, D.: Spherical symmetry as a test case for unconstrained hyperboloidal evolution. Class. Quant. Grav. 32 (17), 175010 (2015) https://doi.org/10.1088/0264-9381/32/17/ 175010 arXiv:1412.3827 [gr-qc]\n- [64] Gasperin, E., Gautam, S., Hilditch, D., nuales, A.: The Hyperboloidal Numerical Evolution of a Good-Bad-Ugly Wave Equation. Class. Quant. Grav. 37 (3), 035006 (2020) https://doi.org/10.  \n1088/1361-6382/ab5f21 arXiv:1909.11749 [gr-qc]  \n- [65] York, J.W. Jr.: Conformal 'thin sandwich' data for the initial-value problem. Phys. Rev. Lett. 82 , 1350-1353 (1999) https://doi.org/10.1103/PhysRevLett.82.1350 arXiv:gr-qc/9810051\n- [66] Pfeiffer, H.P., York, J.W. Jr.: Extrinsic curvature and the Einstein constraints. Phys. Rev. D 67 , 044022 (2003) https://doi.org/10.1103/PhysRevD.67.044022 arXiv:gr-qc/0207095\n- [67] Garfinkle, D., Mead, L.: Cosmological initial data for numerical relativity. Phys. Rev. D 102 (4), 044022 (2020) https://doi.org/10.1103/PhysRevD.102.044022 arXiv:2006.16360 [gr-qc]\n- [68] Bentivegna, E.: Solving the Einstein constraints in periodic spaces with a multigrid approach. Class. Quant. Grav. 31 , 035004 (2014) https://doi.org/10.1088/0264-9381/31/3/035004 arXiv:1305.5576 [gr-qc]\n- [69] Baumgarte, T.W., Murchadha, N.O., Pfeiffer, H.P.: The Einstein constraints: Uniqueness and nonuniqueness in the conformal thin sandwich approach. Phys. Rev. D 75 , 044009 (2007) https://doi. org/10.1103/PhysRevD.75.044009 arXiv:gr-qc/0610120\n- [70] Walsh, D.M.: Non-uniqueness in conformal formulations of the Einstein constraints. Class. Quant. Grav. 24 , 1911-1926 (2007) https://doi.org/10.1088/0264-9381/24/8/002 arXiv:gr-qc/0610129\n- [71] York, J.W. Jr.: Kinematics and Dynamics of General Relativity. In: Workshop on Sources of Gravitational Radiation (1978)\n- [72] Corman, M., East, W.E.: Starting inflation from inhomogeneous initial conditions with momentum. JCAP 10 , 046 (2023) https://doi.org/10.1088/1475-7516/2023/10/046 arXiv:2212.04479 [gr-qc]\n- [73] Aurrekoetxea, J.C., Clough, K., Muia, F.: Oscillon formation during inflationary preheating with general relativity. Phys. Rev. D 108 (2), 023501 (2023) https://doi.org/10.1103/PhysRevD.108. 023501 arXiv:2304.01673 [gr-qc]\n- [74] Jong, E., Aurrekoetxea, J.C., Lim, E.A., Fran¸ca, T.: Spinning primordial black holes formed during a matter-dominated era. JCAP 10 , 067 (2023) https://doi.org/10.1088/1475-7516/2023/10/067 arXiv:2306.11810 [astro-ph.CO]\n- [75] Elley, M., Aurrekoetxea, J.C., Clough, K., Flauger, R., Giannadakis, P., Lim, E.A.: Robustness of inflation to kinetic inhomogeneities (2024) arXiv:2405.03490 [astro-ph.CO]\n- [76] Brady, S.E., Arest'e Sal'o, L., Clough, K., Figueras, P., S., A.P.: Solving the initial conditions problem for modified gravity theories. Phys. Rev. D 108 (10), 104022 (2023) https://doi.org/10. 1103/PhysRevD.108.104022 arXiv:2308.16791 [gr-qc]\n- [77] Bamber, J., Aurrekoetxea, J.C., Clough, K., Ferreira, P.G.: Black hole merger simulations in wave dark matter environments. Phys. Rev. D 107 (2), 024035 (2023) https://doi.org/10.1103/ PhysRevD.107.024035 arXiv:2210.09254 [gr-qc]\n- [78] Aurrekoetxea, J.C., Clough, K., Bamber, J., Ferreira, P.G.: Effect of Wave Dark Matter on Equal Mass Black Hole Mergers. Phys. Rev. Lett. 132 (21), 211401 (2024) https://doi.org/10.1103/ PhysRevLett.132.211401 arXiv:2311.18156 [gr-qc]\n- [79] Giblin, J.T., Mertens, J.B., Starkman, G.D.: Departures from the Friedmann-Lemaitre-RobertstonWalker Cosmological Model in an Inhomogeneous Universe: A Numerical Examination. Phys. Rev. Lett. 116 (25), 251301 (2016) https://doi.org/10.1103/PhysRevLett.116.251301 arXiv:1511.01105 [gr-qc]\n- [80] Joana, C., Clesse, S.: Inhomogeneous preinflation across Hubble scales in full general relativity. Phys. Rev. D 103 (8), 083501 (2021) https://doi.org/10.1103/PhysRevD.103.083501 arXiv:2011.12190 [astro-ph.CO]\n- [81] Goldwirth, D.S., Piran, T.: Initial conditions for inflation. Phys. Rept. 214 , 223-291 (1992) https:", 'https://doi.org/10.1016/0370-2693(82)91219-9': 'of cosmological perturbations through a nonsingular bounce. Phys. Rev. D 88 , 083509 (2013) https://doi.org/10.1103/PhysRevD.88.083509 arXiv:1308.3044 [gr-qc]  \n- [170] Allen, L.E., Wands, D.: Cosmological perturbations through a simple bounce. Phys. Rev. D 70 , 063515 (2004) https://doi.org/10.1103/PhysRevD.70.063515 arXiv:astro-ph/0404441\n- [171] Coleman, S.R.: The Fate of the False Vacuum. 1. Semiclassical Theory. Phys. Rev. D 15 , 2929-2936 (1977) https://doi.org/10.1103/PhysRevD.16.1248 . [Erratum: Phys.Rev.D 16, 1248 (1977)]\n- [172] Coleman, S.R., De Luccia, F.: Gravitational Effects on and of Vacuum Decay. Phys. Rev. D 21 , 3305 (1980) https://doi.org/10.1103/PhysRevD.21.3305\n- [173] Steinhardt, P.J.: Natural inflation. In: Very Early Universe, pp. 251-266 (1983)\n- [174] Vilenkin, A.: The Birth of Inflationary Universes. Phys. Rev. D 27 , 2848 (1983) https://doi.org/ 10.1103/PhysRevD.27.2848\n- [175] Linde, A.D.: ETERNAL CHAOTIC INFLATION. Mod. Phys. Lett. A 1 , 81 (1986) https://doi. org/10.1142/S0217732386000129\n- [176] Guth, A.H.: Eternal inflation and its implications. J. Phys. A 40 , 6811-6826 (2007) https://doi. org/10.1088/1751-8113/40/25/S25 arXiv:hep-th/0702178\n- [177] Aguirre, A.: Eternal Inflation, past and future (2007) arXiv:0712.0571 [hep-th]\n- [178] Aguirre, A., Johnson, M.C.: A Status report on the observability of cosmic bubble collisions. Rept. Prog. Phys. 74 , 074901 (2011) https://doi.org/10.1088/0034-4885/74/7/074901 arXiv:0908.4105 [hep-th]\n- [179] Kleban, M.: Cosmic Bubble Collisions. Class. Quant. Grav. 28 , 204008 (2011) https://doi.org/10. 1088/0264-9381/28/20/204008 arXiv:1107.2593 [astro-ph.CO]\n- [180] Carone, C.D., Guth, A.H.: THE DYNAMICS OF THICK DOMAIN WALLS IN AN INHOMOGENEOUS INFLATIONARY MODEL. Phys. Rev. D 42 , 2446-2457 (1990) https://doi.org/10. 1103/PhysRevD.42.2446\n- [181] Hawking, S.W., Moss, I.G., Stewart, J.M.: Bubble Collisions in the Very Early Universe. Phys. Rev. D 26 , 2681 (1982) https://doi.org/10.1103/PhysRevD.26.2681\n- [182] Blanco-Pillado, J.J., Bucher, M., Ghassemi, S., Glanois, F.: When do colliding bubbles produce an expanding universe? Phys. Rev. D 69 , 103515 (2004) https://doi.org/10.1103/PhysRevD.69.103515 arXiv:hep-th/0306151\n- [183] Takamizu, Y.-i., Maeda, K.-i.: Collision of domain walls in asymptotically anti de Sitter spacetime. Phys. Rev. D 73 , 103508 (2006) https://doi.org/10.1103/PhysRevD.73.103508 arXiv:hepth/0603076\n- [184] Takamizu, Y.-i., Kudoh, H., Maeda, K.-i.: Dynamics of colliding branes and black brane production. Phys. Rev. D 75 , 061304 (2007) https://doi.org/10.1103/PhysRevD.75.061304 arXiv:grqc/0702138\n- [185] Freivogel, B., Horowitz, G.T., Shenker, S.: Colliding with a crunching bubble. JHEP 05 , 090 (2007) https://doi.org/10.1088/1126-6708/2007/05/090 arXiv:hep-th/0703146\n- [186] Takamizu, Y.-i., Chernoff, D.: Collisions of false vacuum bubbles in cylindrical symmetry (2017) arXiv:1710.02701 [gr-qc]\n- [187] Johnson, M.C., Peiris, H.V., Lehner, L.: Determining the outcome of cosmic bubble collisions in full General Relativity. Phys. Rev. D 85 , 083516 (2012) https://doi.org/10.1103/PhysRevD.85.083516 arXiv:1112.4487 [hep-th]  \n023520 arXiv:0709.3510 [astro-ph]', '1361-6382/aba28b arXiv:2002.05177 [gr-qc]': '- [319] Macpherson, H.J.: Cosmological distances with general-relativistic ray tracing: framework and comparison to cosmographic predictions. JCAP 03 , 019 (2023) https://doi.org/10.1088/1475-7516/ 2023/03/019 arXiv:2209.06775 [astro-ph.CO]\n- [320] Grasso, M., Villa, E.: BiGONLight: light propagation with bilocal operators in numerical relativity. Class. Quant. Grav. 39 (1), 015011 (2022) https://doi.org/10.1088/1361-6382/ac35aa arXiv:2107.06306 [gr-qc]\n- [321] Macpherson, H.J., Heinesen, A.: Luminosity distance and anisotropic sky-sampling at low redshifts: Anumerical relativity study. Phys. Rev. D 104 , 023525 (2021) https://doi.org/10.1103/PhysRevD. 104.023525 arXiv:2103.11918 [astro-ph.CO]. [Erratum: Phys.Rev.D 104, 109901 (2021)]\n- [322] Williams, M.J., Macpherson, H.J., Wiltshire, D.L., Stevens, C.: First investigation of void statistics in numerical relativity simulations (2024) arXiv:2403.15134 [astro-ph.CO]\n- [323] Koksbang, S.M., Heinesen, A., Macpherson, H.J.: Redshift drift in a universe with structure III: Numerical relativity (2024) arXiv:2404.06242 [astro-ph.CO]\n- [324] Tian, C., Anselmi, S., Carney, M.F., Giblin, J.T., Mertens, J.B., Starkman, G.: Question of measuring spatial curvature in an inhomogeneous universe. Phys. Rev. D 103 (8), 083513 (2021) https://doi.org/10.1103/PhysRevD.103.083513 arXiv:2010.07274 [astro-ph.CO]\n- [325] Oestreicher, A., Koksbang, S.M.: Backreaction in Numerical Relativity: Averaging on Newtonian gauge-like hypersurfaces in Einstein Toolkit cosmological simulations (2024) arXiv:2408.03049 [astro-ph.CO]'}