{'fraction_non_alphanumeric': 0.029753835176596504, 'fraction_numerical': 0.00927577595433464, 'mean_word_length': 4.624197431781702, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integration of inhomogeneous cosmological spacetimes in the BSSN formalism

We present cosmological-scale numerical simulations of an evolving universe in full general relativity (GR) and introduce a new numerical tool, {\sc CosmoGRaPH}, which employs the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism on a 3-dimensional grid. Using {\sc CosmoGRaPH}, we calculate the effect of an inhomogeneous matter distribution on the evolution of a spacetime. We also present the results of a set of standard stability tests to demonstrate the robustness of our simulations.


I. INTRODUCTION
Current cosmological work typically relies on a perturbative approach or a Newtonian-gravity approximation [1,2]. Although such approximate methods have led to very precise simulations, as we acquire and anticipate sub-percent-level precision measurements of the structure in and of the present-day universe, more accurate simulations and therefore more accurate methods are needed. This is not only to provide more accurate results, but also to yield insights into physical processes not previously appreciated.
In this paper we present a numerical implementation of the fully general-relativistic BSSN formalism for studying gravitational systems [3,4]. This formalism is most often applied in regimes of strong gravity -compact object dynamics. However, it has also been applied to cosmological problems in the early universe, such as critical collapse in a radiation fluid [5,6], black hole lattices (eg. [7]), and more (for a recent review, see [8]). Here we examine the performance of the BSSN formalism in the context of a matter-dominated cosmological spacetime, and draw comparisons between FLRW spacetimes and perturbed FLRW spacetimes.
The equations we evolve are fully nonlinear parameterizations of GR, and therefore formidable to work with analytically -one reason approximations are commonly made. Nevertheless, the nonlinear terms are few enough that, depending on gauge choice, numerically integrating the full unconstrained Einstein equations does not require significantly more computational resources than working in a linearized gravity regime. In this paper we make the case that nonlinear gravitational effects can be taken into account in cosmological simulations using such a formalism, studying whether a fully unconstrained, general relativistic model of a matter-dominated inhomogeneous universe is consistent with FLRW approximations.
Although the system that we study does not depend on an averaging scheme or background solution, there are several assumptions we make in order to simplify numerical integration. The most significant of these are that we use a pressureless fluid with zero initial coordinate velocity, limiting the dynamics we can probe in the matter sector; the simulation contains a short-wavelength cutoff due to finite available resolution; and the system contains a long-wavelength cutoff due to periodic boundary conditions. Here we examine a numerical spacetime beyond the FLRW approximation in synchronous gauge subject to these constraints, obtaining a full nonlinear solution in a cosmological context. To our knowledge, this is the first time that such a formalism has been applied to an unconstrained, inhomogeneous, 3+1-dimensional cosmological spacetime with a pressureless fluid.
In this work we first describe a simple way to determine initial conditions for a universe containing matter and a cosmological constant term. We then numerically evolve a system described by these initial conditions, using a mildly modified version of the standard BSSN equations. We discuss the subsequent evolution of the system, limitations we encounter, and some potential directions for future work.

A. The BSSN Formalism
The BSSN equations are derived from the original ADM formalism [9] of GR, and provide a numerically stable scheme for evolving Einstein's equations. (For a full textbook treatment, see [10].) In this formulation, the metric is written as where α and β i are referred to as the 'lapse' and the 'shift' respectively. The BSSN equations evolve the 3-metric for a particular gauge choice, along with the extrinsic curvature K ij . The metric is rescaled by a conformal factor, γ ij ≡ e 4φγ ij , with det(γ ij ) = 1. The extrinsic curvature is decomposed into a trace, K, and a conformally related trace-free part,Ā ij , whose indices are raised and lowered by the conformal metricγ ij ; arXiv:1511.01106v3 [gr-qc] 24 Jun 2016 The full dynamical system is The source terms written in terms of the stress-energy tensor T µν are where n µ = (−α, 0) and S = γ ij S ij . Our software, the Cosmological General Relativistic and Perfect-Fluid Hydrodynamics (CosmoGRaPH) code, allows for an arbitrary lapse and shift. However, the work we present here largely makes use of synchronous gauge, or geodesic slicing. In this gauge, the lapse is a fixed constant (α = 1) and there is no shift (β i = 0). With this choice, the evolution equations for the metric reduce to (13) For a flat FLRW solution to Einstein's equations, the BSSN variables can be directly translated to the FLRW metric functions. The FLRW spatial metric is γ ij = a 2 δ ij , meaning γ = det γ ij = a 6 . Inserting this relationship into the BSSN equations and applying the gauge choice will give us a translation between BSSN and FLRW parameters: H ∼ − 1 3 K and a ∼ γ 1/6 = e 2φ . These will be useful for computing initial conditions, and for making comparisons with standard results.
In the BSSN formalism, auxiliary conformal connection variables,Γ withΓ ijk =γ ilΓ l jk , are evolved simultaneously with the metric functions. These are used to eliminate terms with mixed derivatives when calculating the Ricci tensor. One writes whereD is the covariant derivative associated withγ ij .
The equation of motion forΓ i is (18) While it might seem odd that we evolve so many (seemingly redundant) parameters, this strategy helps provide numerical stability in the BSSN framework. By evolving the metric coefficients, extrinsic curvature and connections separately, we can check for consistency throughout the simulations. These constraints are referred to as the Hamiltonian and momentum constraints, and should be obeyed on spacetime slices both initially and throughout the simulation. In terms of BSSN variables the Hamiltonian and momentum constraints are respectively and In order to quantify the amount of constraint violation relative to the energy scales involved in the problem, we compare the calculated constraint violation to the root sum of squares of the individual terms in each respective constraint equation -the "scale" of the constraint. We label these two denominators [H] and [M], with M = |M i |.
Alongside the metric, we evolve a pressureless fluid describing the matter content of the universe. Our code was written to utilize a flux-conservative form of the relativistic hydrodynamic equations [11]; however, here we are primarily interested in results for a w = 0 cosmological fluid with rest-mass density ρ 0 with no initial velocity component. The equation of motion for the matter fluid in the absence of any initial velocity is then a simple conservation law, allowing us to optimize the fluid code to model a static source, The contributions of the matter components to the source terms are then ρ = ρ 0 + ρ Λ , S i = 0, S ij = −γ ij ρ Λ , and S = −3ρ Λ .
We fix physical scales by setting our grid side length to L = nH −1 I = N ∆x, where n is an arbitrary fraction and H I is the Hubble scale of an FLRW solution on the initial slice. Combining this with the number of points per grid side, N , allows us to determine the grid spacing ∆x. In programming units, we work in units of the Hubble scale of an analogous FLRW solution at the time the initial slice is set. This also fixes the density of the FLRW solution to be ρ F LRW = (3/8π) (n/N ) 2 . As in almost all simulations, numerical inaccuracies accumulate over time leading to substantial growth of constraint violation at late times. In order to alleviate this problem, we have included a diffusive term in the evolution equations as suggested by [12], with varying degrees of strength. These modifications take the form These additional terms do suppress constraint violation at some level, however late-time growth remains large, ostensibly due to the matter coupling. Although we do not explore this idea thoroughly, we do present some preliminary analysis in Sec. III B. Because of this instability, and because we are interested in percent-level effects, we generally restrict ourselves to timescales on To describe the statistical behavior of fields as they evolve, we commonly compute volume-weighted (conformal) averagesf and volume-weighted standard deviations of various fields (f ) on spatial slices.

II. COSMOLOGICAL INITIAL CONDITIONS
The initial surface from which we evolve should satisfy the Hamiltonian and momentum constraint equations. In the case of an FLRW solution, the Hamiltonian equation is one of the Friedman equations, and all terms in the momentum constraint equation are zero. We are interested in solutions similar to this, but we also wish to examine a universe with non-uniform matter density. We do not attempt to specify a matter distribution that is perfectly analogous to our universe, but instead restrict ourselves to initial conditions that are easier to explore in the context of the BSSN formalism.
The Hamiltonian and momentum constraint equations by themselves do not uniquely specify a metric; more degrees of freedom exist than can be specified by the constraint equations. Additionally, even when using choices that simplify the equations (such as in the conformal transverse-traceless decomposition) the equations can still be difficult to solve. To this end, we relax the requirement that we need to specify an exact matter source for initial conditions, instead specifying the conformal factor φ itself. We specify a slice with constant extrinsic curvature whose value is approximately determined by the average matter density, and fluctuations in density and φ set by a power spectrum. At large scales (small k) we allow the power spectrum to scale as P k ∼ k 1 , and at small scales (large k) as P k ∼ k −3 [13]. Given a peak scale k * and peak amplitude at this scale P * , we use the power spectrum In order to set these initial conditions, we decompose ρ into two pieces, ρ K sourcing the trace of the extrinsic curvature K and ρ ψ sourcing the conformal factor ψ ≡ e φ , so the total density is ρ = ρ K + ρ ψ . We also impose that the matter is at rest, so S i = 0. We use a conformally flat metric and set the trace-free part of the extrinsic curvature to zero, leaving us to solve two simpler equations: The nonlinear ψ-equation is difficult to solve for a fixed matter source ρ ψ , with attempted relaxation and iterative solution methods tending to find the ψ = 0 solution. We instead opt to specify initial conditions by setting ψ directly, with a power spectrum P ψ k = k −4 P k . We set the monopole term to unity, and subsequently set ρ ψ = ∇ 2 ψ/(−2πψ 5 ). A power spectrum realization can be seen in Fig. 1. We note that because we cannot directly specify that ρ ψ has zero mean, the extrinsic curvature K will not necessarily be set by the average of the total densityρ. We generate a Gaussian random field using the matter power spectrum in Eq. 27, specifying a unitary monopole. The field ρ ψ is then calculated using finite differencing from Eq. 28. When the Hamiltonian constraint is subsequently calculated, this means the predominant error should only be due to finite precision. The matter field is constructed as ρ = ψ 5 ρ ψ + ρ K + ρ Λ . Lastly, the trace of the extrinsic curvature is set to K = − 24π(ρ K + ρ Λ ).

A. Code Structure
Here we present details of the code used to evolve the BSSN equations. The engine of the code is written in c++ and employs several c++11 standard features, along with OpenMP [14] for parallelization, FFTW [15] for calculating power spectra and creating initial conditions, and the HDF5 file format [16] for storing samples of data. The simulations presented here are performed on resources available at Kenyon College [17] and the High Performance Computing System at Case Western Reserve, on single nodes that range between 8-64 cores and 64-512 GB of RAM.
The code is divided into different classes each of which evolves a different component. There are two main classes: one that solves the gravitational equations and one that evolves the hydrodynamical equations that define the matter content. There are three additional classes that implement the static w = 0 perfect fluid examined here, a class that defines a cosmological constant, and another that can integrate the pure Friedmann equations for a perfect fluid.
The most important of these is the class that evolves the gravitational equations; this "BSSN class" provides functionality to evolve the BSSN fields in the presence of an arbitrary matter source using a standard fourthorder Runge-Kutta integrator, with four memory registers per gravitational field and an additional register for storing source terms. We informally note that this code runs at a speed comparable to ETK/CACTUS [18,19], but seems to consume roughly an order of magnitude less memory. The BSSN class also includes routines for performing computations such as calculating average momentum and Hamiltonian constraint violation. The results we present are almost exclusively in synchronous gauge, and as such the code has been optimized to take advantage of its simplicity. However, the BSSN class is still capable of evolving the metric in an arbitrary gauge.
The other major class is the "Hydrodynamics class," which allows us to evolve a perfect fluid with a cosmological equation of state. This class has been simplified for a static w = 0 fluid, with further work needed to implement an arbitrary equation of state. As either an N-body code or special lattice methods are required to evolve a pressureless fluid with a non-zero velocity component, we constrain ourselves to work with initial conditions where the fluid is static. In the chosen gauge, the equations of motion simplify such that the conserved fluid density does not evolve, ∂ tD = 0.
Finite-difference stencils are implemented up to eighthorder, as numerical stability and accuracy is highly dependent on structure at smaller scales, as demonstrated below. Unless otherwise specified, the stencils used are generally O(∆x 6 ).
The simulations have periodic boundary conditions, imposed at the level of an array indexing function. However, alternate boundary conditions can be implemented. With our code we expect to be able to run for an arbitrarily large grid-although it is currently limited to singlenode operation. We have been able to run at resolutions up to N 3 = 768 3 on a single 512 GB RAM node.

B. Standard Code Tests: Analytic Cases
We present several code tests where we evaluate several known, analytic solutions, intended to examine the behavior of fields that will dominate the dynamics in cosmological scenarios. The intention is to examine the validity of results presented, and to explore in what regimes the code will produce accurate and interesting results.
We begin with two standard "Apples with Apples" tests [20], the robust stability test and linearized wave test.

Robust Stability Test
The robust stability test examines the growth of small amounts of numerical error by looking at the evolution of noise around a Minkowski background in a simulation with no matter content. Fig. 2 displays this deviation from Minkowski space with different resolutions. In this test, we use an asymmetric box in which there are only 6 points in the y-and z-directions and a variable number of points in the x-direction. In [20], the number of points in the x-direction is parameterized by ρ AwA (the subscript here is added to avoid confusion with the density parameter) such that the number of points in the x-direction is N x = 50ρ AwA .

Linearized Wave Test
The linearized wave test examines the behavior of linear metric fluctuations after a large number of box lightcrossing times. Since we have a plane wave, we work with the asymmetric box of the previous section; where N y = N z = 6 and N x = 50ρ AwA . Spacetime is initialized with plane wave traveling in the x-direction-the direction with the largest resolution-with an amplitude small enough that nonlinear effects are at the level of roundoff error. Fig. 3 shows the difference between the numeri- The matter coupling can be seen to introduce a growing mode; the fixed error present in the staticD field will source the metric continuously. The constraint damping term helps diminish constraint violation short-term. cal result of this wave compared to the expectation from linearized gravity.

Black Hole Profile
To verify the validity of the code for fluctuations in the conformal field φ and gauge fields, we examine an analytic black hole solution with a "trumpet" geometry [21] in the '1 + log' slicing. The solution is stable at remarkably low resolution, N 3 = 32 3 , and with second-order finite-derivative stencils. After initializing our grid with the solution, we see it relax slightly, and subsequently remain stable. This solution is run in vacuum, using a time-independent maximal slicing gauge condition for the lapse, and hyperbolic Gamma driver shift condition [22]. At the outer boundary we naïvely apply fixed boundary conditions. In the interior of the black hole the lapse approaches zero, so the value of the conformal field φ is set such that the derivatives of φ at nearby points will be calculated to a good approximation. We plot the conformal field, φ, in the simulation and compare it to the analytic black hole solution in Fig. 4.

FIG. 4:
A cross section of the conformal field φ after a period of relaxation, analytically (dashed green line) and numerically (blue dots) for a black hole.

FRW universe
We also examine the case of a homogeneous flat FLRW universe. For this we simply choose geodesic slicing, or synchronous gauge, with periodic boundary conditions on a very small (N 3 = 8 3 ) grid. Since there are no spatial inhomogeneities in the simulation, spatial resolution should not play a role in its validity. Indeed, Fig. 5 shows that the software can reproduce FLRW cosmology to high accuracy.

C. Cosmological Code Tests: Variation of Parameters
To test our code in a cosmological regime, we define a fiducial model and then vary physical and numerical parameters to test the stability and reliability of our simulations. The fiducial model has N 3 = 128 3 points on each side; each side of the simulation is half the corresponding FLRW Hubble scale, L = H −1 I /2. We employ O(∆x 6 ) finite-derivative stencils, and ∆t = ∆x/10. The initial conditions have a power spectrum cut at k cutoff = 10/128 ∆k, peak frequency of k ∼ 7/128 ∆k, and peak amplitude such that σ ρ /ρ ∼ 0.04. We note that the resolution we have chosen corresponds to ∆x ∼ 16 MPc, somewhat close to the 11.5 MPc σ 8 scale. However, even without the cut, the amplitude corresponds to σ ρ /ρ ∼ 0.1. Thus a combination of the spectrum cutoff and smaller spectrum amplitude result in relatively small amplitude fluctuations compared to the measured value σ 8 ∼ 0.8 today [23]; thus we are studying a smoother universe.
So that we have a method of comparison, we present some baseline information about our fiducial model. Fig. 6 shows how well our fiducial model satisfies the physical constraints, Eqs. 19 and 20, and Fig. 7 shows how inhomogeneities of the extrinsic curvature, σ K /K, are generated over the course of the simulation.

Constraint Violation
Here we explore the effect of the constraint damping terms from Eq. 22. Increasing the magnitude of c H helps reduce the maximum amount of Hamiltonian constraint violation, although it does not have much effect on the overall growing mode. Fig. 8 shows the effect of this term on level of Hamiltonian constraint violation and on the generation of inhomogeneities of the extrinsic curvature,

Varying the Resolution
Here we examine the effects of changing the resolution while keeping the same initial conditions (picking initial conditions with identical amplitudes for k). We see that the amount of constraint violation diminishes as the resolution is increased, something expected due to the increasing smoothness of the fields, making derivatives more accurate. Fig. 9 shows that we see no qualitative difference in the generation of inhomogeneities and marginally better satisfaction of the Hamiltonian constraint, Eq. 19, as we increase the resolution of the simulations.
A more precise test of resolution variation can be seen in Fig. 10. We should be able to predict how much better the constraints can be satisfied as we increase the resolution; specifically, as we double the number of points on a side, we should see the ratio of errors be O(∆x 6 )/O(∆x 6 /2 6 ) ∼ 2 6 .
We have implemented a number of finite-differencing stencils, seeing reduced error for higher-order stencils as expected. The use of O(dx 2 ) stencils in particular lead to unacceptable growth of constraint violation in our simulations. We also tested changing the time resolution ∆t to check for convergence. The results presented here agree with runs for ∆t an order of magnitude smaller.

Varying the Spectrum Cutoff
As a final technical test, we present results as the spectrum cutoff is varied. As the cutoff approaches smaller | k|, we see that the amplitude of fluctuations is significantly diminished, as expected, when modes are removed. The relative amount of constraint violation also diminishes as the fields become increasingly smooth. As the cutoff is raised, we see the opposite happen: structure at increasingly small (and unresolvable) scales is created, and constraint violation increases. The amount of constraint violation after an e-fold is a part in ten when c = 40, and is of order a part in 10 11 for c = 2.
This is an important test, specifically when compared to the resolution test, Sec. III C 2, as it shows that we have control over the finite differencing scheme; for example, doubling the value of k cutoff resolves the same physical modes as would halving the physical size of the box (and keeping the resolution constant), or doubling the resolution of the box. While the numerical strategy differs for these three cases by many factors of ∆x, ∆t, N , etc. we see excellent agreement in the results from our code. We have demonstrated that these three situations produce the same physical results and show that we are using our software well within the regime of its validity.

IV. RESULTS
Now that we have demonstrated the stability and reliability of our simulations, we complement the main results of these simulations as presented in [24]. We are primarily interested in understanding how inhomogeneous matter generates local differences in the expansion rate. We do this in two steps. First, we look to see if we generate fluctuations of the extrinsic curvature, σ K /K, which, in our gauge, is is the source for the conformal factor, φ. Once this is established, we will run two tests to tell whether or not this generation is an artifact of our scheme or has a measure that shows local deviation from FLRW.
We compare the fiducial model to the corresponding FLRW simulation with the same initial average initial density,ρ. The results are shown in Figs. 12 and show that average quantities are very consistent with those in the corresponding homogeneous FLRW universe. Average parameters correspond very well with FLRW parameters for the w = 0 case we study here. However,  2,4,8,10,12,14,16,20,40, from blue to red (bottom to top). Note that as we initialize more modes, we are also increasing σρ/ρ.
we do begin to see small deviations from the FLRW approximation that are robust to resolution increase. The significance of this result is on the order of the constraint violation.
On the other hand, we always generate deviations from homogeneity in the extrinsic curvature. We already presented this result for our fiducial model in Fig. 7. As the simulation progresses we see significant local departures from constant K.
Of course, we want to explore the relationship between the growth of inhomogeneities of K with increasingly inhomogeneous matter sources. Varying the amplitude has a corresponding effect on the matter over-density (analogue of σ 8 ).
Now that we have established the existence of spatial variations of the extrinsic curvature, K, we seek to explore the behavior of additonal measures associated with the spacetime. The first test is to see if there is a departure from FLRW predictions for the proper lengths of paths; in pure FLRW cosmology, the proper length of any path will depend on its initial length and the scale factor, a FLRW . Since we use synchronous gauge, it is a good test to chose arbitrary paths and calculate the ratio of the lengths of these paths as a function of time-if they vary from the FLRW prediction, there is a clear departure from the homogeneous model. We calculate the length of paths by between two arbitrary points, x 1 and x 2 . Fig. 14 and has solutions k µ . In a flat FLRW universe, the Killing vectors associated with translations are (k µ x , k µ y , k µ z ) = ((0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)). In an inhomogeneous universe, the killing equation should not be perfectly satisfied. We look at the contracted killing equation associated with translations in the i-direction, We look at the growth of this quantity relative to the timescale on which the metric evolves, K.

V. DISCUSSION
We have shown that non-linear GR effects are important on sub-Hubble scales in the presence of inhomogeneities in the universe. Without assuming anything about the dynamical system (e.g. not using symmetry to simplify the dynamical GR equations) we have performed a simulation of the universe that does not rely on a background or averaging procedure. Using this method we generate FLRW behavior for homogeneous cosmologies and see deviation from this for inhomogeneous ones. We show that, even beginning with constant extrinsic curvature, K, inhomogeneities will generate local fluctuations in this parameter. The color coding goes from paths of initial coordinate distance of 4∆x (purple) to paths of initial coordinate distance 128∆x (red) in spectral order. The bottom panel shows the same set of paths, but compares the ratio of proper length to initial proper length with the FLRW estimate.
Physical observables need to be addressed independently. For example, the magnitude of the effect on a Hubble diagram would require the integration of photon geodesics during the simulation. We defer this, and other observable tests, for future work.
As the dominant contribution to the BSSN equations in an FLRW spacetime is simply the FLRW solution, it is in principle possible to subtract the FLRW solution from the BSSN equations, and evolve variables representing differences between the FLRW solution and an inhomogeneous solution. Note that this does not constitute an approximation: the equations are still fully relativistic under this procedure. We do not explore this idea here, however we may in future work.

VI. ACKNOWLEDGMENTS
We would like to express a special thanks to Thomas Baumgarte for extremely valuable conversations and advice during early stages of this project. JTG would like to thank Eugene Lim and Katy Clough for useful discus-sions and for the hospitality at King College London at which these conversations took place. JBM would like to thank Shaden Smith for initial contributions to the code. JBM and GDS are supported by a Department of Energy grant DE-SC0009946 to the particle astrophysics theory group at CWRU. JTG is supported by the National Science Foundation, PHY-1414479. We acknowledge the National Science Foundation, the Research Cor-