Signatures of primordial gravitational waves in matter power spectrum

We simulate the evolution of a dust universe from $z=1089$ to $z=0$ by numerically integrating the Einstein's equation for a spatially flat Friedmann-Lemaire-Robertson-Walker (FLRW) background spacetime with scalar perturbations which are derived from the matter power spectrum produced with the Code for Anisotropies in the Microwave Background (CAMB). To investigate the effects of primordial gravitational waves (GWs) on the inhomogeneity of the universe, we add an additional decaying, divergenceless and traceless primordial tensor perturbation with its initial amplitude being $3\times 10^{-4}$ to the above metric. We find that this primordial tensor perturbation suppresses the matter power spectrum by about $0.01\%$ at $z=0$ for modes with wave number similar to its. This suppression may be a possible probe of a GWs background in the future.


A. Initial conditions for background
Since we will perform large-scale cosmological simulations instead of the simulations of black-hole-binary-like astrophysical system, we modify the file EOS Omni Module.F90 in Einstein Toolkit to replace the default unit system: M = G = c = 1 with the new one: 1Mpc = G = c = 1 [33]. Under this new unit system and with the cosmological parameters consistent with Planck 2018 results [34] as shown in Tab. I, the matter density of our universe is ρ P m = 6.0 × 10 −9 × 0.3166 at z = 0, henceρ P m = 6.0 × 10 −6 × 0.3166 at z = 9. Considering a fiducial universe whose matter densityρ F m (z) is equal toρ P m (z), the scale factor of this fiducial one a F = 10a P as shown in Tab. II means that the comoving matter density of it isρ F m * = 6.0 × 10 −6 × 0.3166 as shown in Tab. III. Here we will turn to a blown-up fiducial universe by 109 times to mimic our universe in simulations: setting the scale factor used during simulations as a S = 109a F as shown in Tab. II and keeping the comoving matter density beingρ S m * = 6.0 × 10 −6 × 0.3166 as shown in Tab. III. That is to say, simulations withρ S m * = 6.0 × 10 −6 × 0.3166 from a S = 1 to a S = 1090 can give the evolution of our universe withρ P m * = 6.0 × 10 −9 × 0.3166 from a P = 0.00092 to a P = 1 when we analyze the results from simulations taking this blowing-up by 109 times into consideration and regardless of the existence of dark energy and radiation.   Three conversions between scale factor and redshift z. a P follows the usual convention in cosmology. a S is used during our simulations. a F is a fiducial one which relates the former two.
All in all, we set the initial scale factor and matter background density for simulations as a S init = 1 andρ S m,init = 6.0 × 10 −6 × 0.3166 respectively.

B. Initial conditions for perturbations
In the conformal Newtonian gauge, the line element that includes both the scalar and tensor perturbations to a spatially flat FLRW background spacetime is where η is the conformal time, δ ij is the identity matrix, Ψ is the Newtonian potential, Φ the spatial curvature perturbation and h ij is a divergenceless, traceless and symmetric tensor. At the beginning of simulations, it's reasonable to take (1) as the universe's metric and transfer it into the form of (3 + 1) formalism where α is the lapse function which satisfies the harmonic slicing here: ∂ t α = − 1 3 α 2 K, β i is the shift vector which is set as β i = 0 here and γ ij is the spatial metric which evolves depending on the extrinsic curvature K ij as (∂ t −L β )γ ij = −2αK ij . Therefore, the initial data for thorn ADMBase and HydroBase can be derived from the solutions at η = 0 to Einstein's equation for (1).
Given the energy-momentum tensor of a perfect fluid without the anisotropic stress tensor T µν = (ρ+P )u µ u ν +P g µν , we can give the evolutions of a S andρ S m according to the dust (P ρ ≡ρ S m (1+δ)) solutions to the zero-order Einstein equations for (1) It's obviously that a S ,ρ S m , ξ and η are functions of t for FLRW background spacetime and they will become spacedependent in an inhomogeneous spacetime. For the latter case, we still take them as background quantities by taking the average of them across the simulation box. Also, we can give the evolutions of perturbations according to the solutions to first-order Einstein equations for (1) where ε s ij with s = ×, + are transverse and traceless polarization tensors and each of h s k (η) evolves independently and satisfies h s . According to (3) and (4), at η = 0 (or ξ = 1), the initial data will dependent on a S init ,ρ S m,init , f (x i ), h s k (0) and η 0 . The last plot in Fig. 1 shows the distribution of spatial curvature perturbations Φ(x i ) (or f (x i )) at a S = 1. In fact, we use the function make gaussian random field in c2raytools [35] to generate the density perturbations δ(x i ) (the second plot in Fig. 1) from the matter power spectrum at z + 1 = 1090 (the first plot in Fig. 1) produced by the Code for Anisotropies in the Microwave Background (CAMB) [36] with parameters listed in Tab. I first. And then we derive f (x i ) from δ(x i ) according to the Fourier version of (4), hence Φ(x i ), Ψ(x i ) and v i (x i ) (the third plot in Fig. 1). As for tensor perturbations, we here only consider one single mode with k = 2π L and the space distribution as cos( 2π L z), where L = 1000 is the length of one side of our simulation box with x i in [−500, 500]. And we set its initial amplitude h s 2π L (0) = 10 −3 , but it has crossed inside the horizon and decayed by 70% when η 0 2

III. RESULTS
Our simulations are performed at 160 3 , 80 3 and 40 3 resolution and end at a S = 1090. Due to the coincidence of the black curve drawn by 3j1[k(η0+η)] is the spherical Bessel functions of order one) and the red one which is the evolution of γ12(η) (a S ) 2 given by simulations with only tensor perturbations, in Fig. 2, we relate γ12(η) (a S ) 2 to the evolution of tensor perturbation h × (η 0 + η) in our simulations. Although, as shown in Fig. 2, there are slight deviations between the red curve and the green one which is the evolution of γ12(η) (a S ) 2 given by simulations with scalar and tensor perturbations, we keep this relation standing. For probing the effects of primordial tensor perturbations on the inhomogeneity of the universe, it's naive to compare the distribution of δ(x i ) at a S = 1090 given by simulations with scalar and tensor perturbations and their counterparts with only scalar perturbations directly, as shown in Fig. 3. Here we will turn to the the matter power spectrum, which is an important statistical quantity and can be detected by many experiments [17][18][19][20]. In the left plot of Fig. 4, the red, blue and cyan curves are the matter P(k) There are four curves in the first plot: the black one is the matter power spectrum at z + 1 = 1090 produced by CAMB [36] with parameters listed in Tab. I; the red, blue and cyan curves are the matter power spectra drawn from density perturbations δ(x i ) by the function power spectrum 1d in c2raytools[35] at 160 3 , 80 3 and 40 3 resolution respectively, while δ(x i ) is generated by the function make gaussian random field in c2raytools from the black curve. The second plot is the distribution of δ(x i ) at a S = 1 and 160 3 resolution. The last two plots are the distribution of v x (x i ) and Φ(x i ) at a S = 1 respectively, which are derived from δ(x i ) according to the Fourier version of (4).
power spectra drawn from density perturbations δ(x i ) at a S = 1090 by the function power spectrum 1d in c2raytools at 160 3 , 80 3 and 40 3 resolution respectively, where δ(x i ) is given by simulations with only scalar perturbations. When taking the tensor perturbations into consideration, we can get similar matter power spectra. Comparing them with the formers, we can find an obvious suppression of matter power spectra for modes with wave number similar to the tensor perturbations', as shown in the right plot of Fig. 4. And comparing the suppression at 160 3 , 80 3 and 40 3 resolution, we can find this suppression converge to about 0.01% if the initial amplitude of the tensor perturbations is 3 × 10 −4 . Even though the initial conditions derived from the matter power spectrum at z + 1 = 1090 satisfy the perturbed Einstein equations, it's still necessary to check that to what extend do these initial data satisfy the Hamiltonian at the origin of our simulation box. The black curve is drawn by is the spherical Bessel functions of order one. The red curve is the evolution of γ 12 (η) (a S ) 2 in simulations with only tensor perturbations. The green curve is the evolution of γ 12 (η) (a S ) 2 in simulations with scalar and tensor perturbations. We can see that the black curve and the red one are almost coincide and there are slight deviations between the red curve and the green one. That is to say, we can relate γ 12 (η) (a S ) 2 to the evolution of tensor perturbation h × (η0 + η) in our simulations. It's worth pointing out that although, at the beginning of simulations, our initial data is derived from matter power spectrum at z + 1 = 1090, we ignore the dark energy in our simulations at the late time. So the black curve has a different trend with color ones in the left plot. The red, blue and cyan curves in the right plot explicitly show the suppression of matter power spectra for modes with wave number similar to the tensor perturbations' at 160 3 , 80 3 and 40 3 resolution respectively. And this suppression converge to about 0.01%.
constraint and the momentum constraint. 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 Fig . 5 shows the evolution of L 2 norms of the Hamiltonian constraint violation and the x-component of momentum constraint violation at 160 3 , 80 3 and 40 3 resolution. We can see that the higher resolution, the larger constraint violation. The reason for this abnormal behaviour is that the initial δ(x i ) generated by the function make gaussian random field in c2raytools from the matter power spectrum at z + 1 = 1090 produced by CAMB is resolution-dependent: the higher resolution leads to δ(x i ) with larger wave number; the scalar perturbations on smaller scales have larger amplitude. As pointed out in [33], one can present the convergence of constraint violation explicitly by transferring raw constraint violation to relative one.

IV. SUMMARY AND DISCUSSION
We simulate a dust universe from a S = 1 (or z = 1089) to a S = 1090 (or z = 0) by numerically integrating the Einstein's equation whose solution at a S = 1 is a spatially flat FLRW metric with scalar perturbations which are derived from the matter power spectrum produced with CAMB. Then we add an additional decaying, divergenceless and traceless primordial tensor perturbation with its initial amplitude being 3 × 10 −4 to the metric as shown in Fig. 2. Simulations at 160 3 , 80 3 and 40 3 resolution converge and show that this primordial tensor perturbation suppresses the matter power spectrum by about 0.01% at z = 0 for modes with wave number k ∼ 0.05 as shown in Fig. 4. In the linear perturbation theory, scalar and tensor perturbations are supposed to be totally decoupled. Therefore, this suppression results from the fully relativistic treatment for Einstein equations. Although there are nonlinear structures formed at the end of simulations (a S = 1090) as shown in Fig. 3, their scales are smaller than tensor perturbations'. So this suppression sown before the tensor perturbations died out and amplified with time is still in linear regime. This suppression may be a possible probe of a GWs background in the future only if the matter power spectrum is measured in high enough precision. Undoubtedly, by the time LSST is in full operation, the required precision for detection of such suppression is still far beyond reach. However, this suppression is an unique signature put by primordial GWs.