Numerical relativity investigation of the effects of gravitational waves on the inhomogeneity of the universe

We numerically integrate the Einstein’s equations for a spatially flat Friedmann–Lemaire–Robertson–Walker (FLRW) background spacetime with a spatial curvature perturbation and evolving primordial tensor perturbations using the Einstein Toolkit. We find that although the primordial tensor perturbation does not play an important role in the evolution of the overdensity produced by the scalar perturbation, there is an obvious imprint left by the primordial tensor perturbation on the distribution of the fractional density perturbation in the nonlinear region. This imprint may be a possible probe of a gravitational waves background in the future.


Introduction
Inflation [1,2] predicts that there is a stochastic gravitational waves (GW) background. Therefore, it's possible to test inflation scenario experimentally through detection of such a GW background. So far, the B-mode polarization of the cosmic microwave background (CMB) is the most promising probe of this GW background [3][4][5]. In the future, 21cm HI emission from the dark ages will be a complementary and even more sensitive probe of this GW background [6,7]. Furthermore, there are some not very competitive probes of this GW background, including weak lensing shear [8,9] and other large-scale structure observables [10,11]. The goal of this paper is to study the effect of gravitational waves on the inhomogeneity of the Universe with numerical relativity, thereby proposing a possible probe of a GW background.
Although cosmological principle points out that the Universe is homogeneous and isotropic on large scales and can be described by a Friedmann-Lemaitre-Robertson-Walker (FLRW) model, it is inhomogeneous and anisotropic on scales smaller than ∼ 80 h −1 Mpc [12,13] today. These inhomogeneities can induce nonlinear general relativistic effects a e-mail: wangke@itp.ac.cn which may be detected by forthcoming cosmological surveys [14][15][16]. Moreover, these nonlinear general relativistic effects on small scales may be accompanied by unexpected nonperturbative behavior on larger scales. This "backreaction" is believed to be the reason of the recent cosmic expansion by [17][18][19][20][21][22][23][24]. We will, however, use a restricted simulation setup by using a fixed FLRW background and subjecting perturbations to periodic boundary conditions, which suppresses the large-scale backreaction effect on the global evolution of the background spacetime [25].
Usually, the linear perturbation theory of General Relativity (GR) is used on large scales and Newtonian N-body simulations (or Newtonian gravity) provide a very good approximation on small scales. To study the nonlinear general relativistic effects, one can do general relativistic N-body simulations as [26]. Undoubtedly, the direct numerical integration of Einstein's equation is the only way without any systematic errors and approximation to study the Universe on all scales. The first cosmological work that is fully non-linear, fully relativistic and does not impose symmetries or dimensionalreductions has been done by [27,28] using CosmoGRaPH. For cosmological purpose, soon afterward, [29] turned to the wide-used Einstein Toolkit [30] to integrate Einstein's equation. Macpherson et al. [31] also studied the inhomogeneous cosmology with Einstein Toolkit by developing a new thorn, FLRWSolver.
Here, our work is based on: a new thorn CFLRW-Solver (a C language counterpart of FLRWSolver) which takes the tensor perturbation into consideration and initializes an almost FLRW Universe with scalar and tensor perturbations; a self-developing thorn CFLRWAnalysis which calculates several derived variables including the comoving time, scale factor, the tensor component of perturbation and the distribution of overdensity at the end of each evolution step; the thorn McLachlan [32][33][34] which evolves spacetime using the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism [35][36][37] and the thorn GRHydro which evolves the hydrodynamical system [38][39][40].
This paper is organized as follows. In Sect. 2, we give the evolution equations of the FLRW background spacetime and scalar and tensor perturbations through solving the zeroorder and first-order of Einstein equations and solutions to them. In Sect. 3, we give the initial conditions of the system with small perturbations needed by the thorn CFLR-WSolver. In Sect. 4, we analyse the results of simulations provided by the thorn CFLRWAnalysis. At last, a brief summary and discussion are included in Sect. 5.
In this paper, we adopt the following conventions: Greek indices run in {0, 1, 2, 3}, Latin indices run in {1, 2, 3} and repeated indices implies summation and we are in a geometric unit system with G = c = 1.

Cosmological perturbations
For a spatially flat FLRW background spacetime, the line element is where η is the conformal time, a is the scale factor and δ i j is the identity matrix. In the conformal Newtonian gauge, the line element that includes both the scalar and tensor perturbations to the metric is where is the Newtonian potential, the spatial curvature perturbation and h i j is a divergenceless, traceless and symmetric tensor. And for a perfect fluid without the anisotropic stress tensor, its energy-momentum tensor with density ρ = ρ 0 + ρ 1 , isotropic pressure P = P 0 + P 1 and 4-velocity The Einstein equations relate the spacetime curvature to the energy-momentum tensor as The zero-order Einstein equations give the Friedmann constraint and evolution equation for the FLRW background spacetime where a prime represents a derivative with respect to the conformal time. According to [31], the dust (P ρ) solution to (5) is where a init and ρ 0,init are the values of a and ρ 0 at η = 0 respectively, ξ is the scaled conformal time and ρ 0 * = ρ 0 a 3 is the conserved comoving density. From the first-order perturbed Einstein equations, we derive equations describing scalar metric perturbations as [31,41] [31] also gives the dust (P ρ) solution to (7) for the growing mode as where f (x i ) is an arbitrary function of space, C 1 = a init 4πρ 0 * and C 2 = − a init 6πρ 0 * . From the spatial part of the first-order perturbed Einstein equations, we have a wave equation One can expand the tensor perturbation in plane waves where ε s i j with s = ×, + are transverse and traceless polarization tensors and each of h s k (η) evolves independently and satisfies According to [42], for modes inside the horizon during matter dominated era, the exact solution is where η + η 0 is the comoving size of horizon.

Initial conditions
To integrate Einstein equations, Einstein Toolkit turns to the metric in the form of (3 + 1) formalism where α is the lapse function, β i is the shift vector and γ i j is the spatial metric and evolves depending on the extrinsic curvature K i j as Now we will use CFLRWSolver to initialize an almost FLRW Universe with small perturbations. First, we should relate the groups of gird function for basic spacetime variables in the thorn ADMBase to the variables in (2): where we have set dt = a √ 1+2 α dη and β i = 0. Here we choose the harmonic slicing which describes the evolution of α. Since K is negative in our simulation, this choice of foliation will be with high computational efficiency. Then, we relate the basic variables and grid functions for hydrodynamics evolutions in the thorn HydroBase to the variables in (3): According to (6), (8) and (12), we set the initial conditions for a dust system with periodic boundary conditions at t = η = 0 as where l = 500 is the half length of one side of our simulation box with x i in [−375, 625], 0 is the amplitude of spatial curvature perturbation and h s 2π L (0) is the amplitude of monochromatic primordial tensor perturbation with wave number k = 2π L before horizon-crossing. For simplicity, we set the initial scale factor a init = 1. And we set L < 4000 and ρ 0,init = 10 −6 so that 2πη 0 L > 1 which implies the tensor perturbation has crossed inside the horizon η 0 2 3 8πρ 0,init at the beginning of simulation. To keep the linear approximation remains valid and save the computational time as much as possible, we set 0 = 10 −4 which means the density contrast is about 10 −3 . Here we will run three simulations at resolution 100 3 with h s (0) = 10 −3 to study the effect of tensor perturbation on the inhomogeneity of universe.

Results
We can analyse the effect of tensor perturbation by comparing the outputs of several variables in the CFLRW-Analysis derived from the basic ones in the ADMBase and HydroBase from those simulations performed above directly. Firstly, the relation dt = a √ 1+2 α dη gives from which we can use (6) to obtain the evolution of a and ρ 0 with the comoving time η , hence the evolution of h × (η 0 +η) and ρ 1 ρ 0 (η) where γ 12 (t) and (ρ 1 + ρ 0 )(t) are basic variables in ADM-Base and HydroBase. Since the background quantities including a, ρ 0 and η are slightly space-dependent in an inhomogeneous spacetime, especially in the nonlinear region in a dynamical gauge like (16), here these quantities are the average ones taken across the simulation box, which will introduce an element of inaccuracy and gauge dependence. (0) = 10 −3 (black solid curve). We can see that the scalar perturbations can produce the tensor perturbations due to the nonlinear effects (red solid curve) as pointed by [43], and the evolution of primordial tensor perturbation (black solid curve) follows the exact solution j 1 [k(η 0 +η)] k(η 0 +η) (green solid curve), where j 1 (z) = sin z−z cos z z 2 is the spherical Bessel functions of order one. Figure 2 shows the evolution of ρ 1 ρ 0 at the point (375, 375, 375) of our simulation box without primordial tensor perturbations (left), and the primordial tensor perturbation's contribution to ρ 1 ρ 0 at the point (375, 375, 375) (right). We can see that the simulation result (black solid curve) deviates from the linear analytic solution (red solid curve) due to the nonlinear effects. Moreover, even though the primordial tensor perturbation die off quickly after the horizon-crossing as shown in Fig. 1, their contribution to ρ 1 ρ 0 grow up quickly in nonlinear regions. Figure 3 shows the distribution of ρ 1 ρ 0 on the x-y plane of z = 375 at the beginning (η = 0) and end (η = 12000) of simulation without primordial tensor perturbations. We can see that the locations of the maximum of ρ 1 ρ 0 are fixed at points (−125, −125, 375), (375, −125, 375), (375, 375, 375) and (−125, 375, 375). That is to say there is almost no interaction between the perturbation peaks. The left one of Fig. 4 shows the contribution of primordial tensor perturbation with 1 Since we set the initial conditions of ρ1 ρ0 and h × as − sin( 2π x i 500 )-like and cos[ 2π(z+125) 500 ]-like respectively, here we choose the particular locations for extracting the results where the maximums are obtained. just h × component to ρ 1 ρ 0 , while the right one shows the contribution from primordial tensor perturbation with only h + at the end of our simulation. We can see that the contributions of both cases are too small to modify the right one of Fig. 3. However, primordial tensor perturbation with different component do leave a characteristic imprint on the distribution of ρ 1 ρ 0 : the tensor perturbation with h × enhances the overdensity; the tensor perturbation with h + enhances the overdensity alone the lines (−125, y, 375) and (375, y, 375) but suppresses the overdensity alone the lines (x, −125, 375) and (x, 375, 375).

Summary and discussion
We have performed three simulations using Einstein Toolkit in this paper: The first one gives the evolution of overdensity ρ 1 ρ 0 in a spatially flat FLRW background spacetime with a spatial curvature perturbation; In the next two, we added an evolving primordial tensor perturbation with just h × or h + component to the spacetime and find that these two components leave a characteristic imprint on the distribution of the fractional density perturbation in the nonlinear region. More precisely: The primordial tensor modes do decay rapidly and never leave the linear regime as shown in Fig. 1. Before they die out, however, they have modified the profile of overdensity slightly, and the modification is amplified with time as shown in the right one of Fig. 2. The Fig. 4 shows the final modification in the nonlinear region, so the FIG.4 can only be produced by the nonlinear numerical relativity code which is a necessary here. These imprints may be a possible probe of a GW background in the future. We do try to suggest exactly how or how much these tenor modes may affect observables including 2D angular power spectrum of large-scale structure tracers (as Fabian Schmidt et al. proposed [10,11], but unfortunately here we just study the scalar modes with wave length l = 500 and there is no suitable primordial power spectrum on hand. And it's worth pointing out that much of the infrastructure for our work was developed and the development of tensor perturbations from scalar perturbations was already demonstrated in the MacPherson et al. paper [31]. For simplicity and correctness, our work takes advantage of their configuration for simulation directly. That is to say their work serves as an important tool for our work. On the other hand our work also goes beyond their work: (a) additional primordial tensor modes have been added to the initial conditions and their evolution has been shown in Fig. 1. (b) Their contribution to the evolution of overdensity has been shown in the right one of Fig. 2. (c) Their final imprint on the distribution of the fractional density perturbation in the nonlinear region also has been given in Fig. 4. Since the half wave length of the perturbations used in our paper is 250, we perform the three-point convergence test will converge to 2 n and O will converge at nth order. Figure 5 shows that the c of ρ 1 + ρ 0 remains close to the second-order value 4. In our paper, we give the initial conditions by solving the perturbed Einstein equations with 0 = 10 −4 . So there is a question that whether these initial data satisfy the Hamiltonian constraint and the momentum constraint or not. Given the 3-Riemann scalar (3) R, the covariant derivative associated with the 3-metric D j , and the matter energy and momentum density as measured by the Eulerian observer E and p i , we can specify the form of the Hamiltonian constraint violation and the momentum constraint violation as and (A3) Figure 6 shows the evolution of the maximum of the Hamiltonian constraint violation and the x-component of momentum constraint violation under three resolutions: ( 1000 6.25 ) 3 , ( 1000 12.5 ) 3 and ( 1000 25 ) 3 . We can see that the constraints remain bounded during the evolution and converge at first order with increasing resolution. In particular the strong initial increases and latter spikes shown in Fig. 6 become gentle with increasing resolution.