Non-linear evolution of cosmologies I: methodology
Abstract
We introduce the method and the implementation of a cosmological simulation of a class of metric-variation models that accelerate the cosmological expansion without a cosmological constant and evade solar-system bounds of small-field deviations to general relativity. Such simulations are shown to reduce to solving a non-linear Poisson equation for the scalar degree of freedom introduced by the modifications. We detail the method to efficiently solve the non-linear Poisson equation by using a Newton-Gauss-Seidel relaxation scheme coupled with multigrid method to accelerate the convergence. The simulations are shown to satisfy tests comparing the simulated outcome to analytical solutions for simple situations, and the dynamics of the simulations are tested with orbital and Zeldovich collapse tests. Finally, we present several static and dynamical simulations using realistic cosmological parameters to highlight the differences between standard physics and physics. In general, we find that the modifications result in stronger gravitational attraction that enhances the dark matter power spectrum by for large but observationally allowed modifications. More detailed study of the non-linear effects on the power spectrum are presented in a companion paper.
I Introduction
Explaining the accelerated cosmological expansion has become one of the greatest challenges in modern day physics (see e.g., Frieman et al., 2008). Of the many potential explanations, two classes are popular today: ones that propose new and unseen form of energy with negative pressure, and ones that modify the Einstein-Hilbert action. Of particular interest in recent literature are a class of modifications to the Einstein action by an addition of non-linear functions of the Ricci scalar, , which have been demonstrated to cause accelerated expansion for a range of functions (Carroll et al., 2004; Nojiri and Odintsov, 2003; Capozziello et al., 2003).
In addition to explaining the cosmic acceleration, such modifications must satisfy stringent constraints imposed by solar system tests (Williams et al., 1996; Anderson et al., 1998; Bertotti et al., 2003; Williams et al., 2004) as well as cosmological constraints imposed by recent dark energy measurements (Knop et al., 2003; Riess et al., 2004; Spergel et al., 2007). Recently, a new class of models was proposed by Hu & Sawicki (Hu and Sawicki, 2007a, hereafter HS model) that has the flexibility to reproduce the observed cosmological acceleration as well as to satisfy solar system bounds on deviations from the Newtonian limits of general relativity. In the proposed model, the scalar degree of freedom introduced by the modification becomes massive in regions of high curvature and is hidden by the so-called chameleon mechanism (Khoury and Weltman, 2004a, b).
In such models that evade both the cosmological background test and solar system tests, the next strongest constraint comes from the effects of the modifications to the intermediate lengths scales, in the regime of cosmological large scale structures. Outside of dense dark matter halos, the light scalar degree of freedom produces a long range fifth force that enhances gravitational attraction and thereby potentially change clustering of matter. Such effects were studied in the linear regime (Song et al., 2007a; Hu and Sawicki, 2007a; Pogosian and Silvestri, 2008), in which large deviations of dark matter power spectrum were obtained even for small modifications. Furthermore, an attempt to model some of the non-linear chameleon physics using a parametrized framework on top of a halo model are presented in Hu and Sawicki (2007b). However, those studies remain inconclusive because the effects of non-linear growth of structure and the chameleon mechanism (which requires non-linear growth) are not included in the analyses. For such non-linear study, a realistic N-body simulation of structure formation in cosmologies is needed.
In this paper, we introduce a formalism and an implementation of a cosmological N-body simulation with fully non-linear treatment of modifications. Although somewhat limited in the attainable dynamic range, we show that our implementation can capture the essential non-linear features of the HS model that have large impact on the cosmological large scale structure formation. In a companion paper, we apply our methodology to realistic cosmological simulations and discuss the implications of the model on the large scale structure of the Universe.
This paper is organized as follows. In §II, we review the general properties of cosmology and give details on the particular model of that we choose. The numerical algorithms to solve for the cosmological large scale structure formation are developed in §III, and are tested in §IV. Several test cases under the HS model, highlighting the different aspects of the model, are presented in §V. In addition, a brief discussion of full cosmological simulations with modifications are presented in §V; a more complete discussion will be presented in a companion paper. Finally, conclusions and discussions of possible future extension of our formalism, are given in §VI.
Ii cosmology
The so-called theories are modifications to general relativity (GR) given by the following action:
(1) |
where is the scalar curvature (i.e., the Ricci scalar), is an arbitrary function of , , and is the Lagrangian of the ordinary matter. Setting recovers ordinary GR without a cosmological constant. We follow the standard procedure of setting and throughout the paper unless otherwise noted. The variation of the action with respect to the metric results in the field equation for ,
(2) |
where , is the unmodified Einstein tensor, and is the energy-momentum tensor. We work with the standard assumption that the Universe is filled with a perfect fluid of the form
(3) |
where and are the fluid energy density and pressure and is the four velocity of the fluid. Because we are only concerned with modeling the effects of late time accelerated expansion, we neglect the relativistic energy component and thus set . For the background Friedmann-Robertson-Walker metric, the modified Einstein equation becomes the modified Friedmann equation
(4) |
where is the Hubble parameter, is the second derivative of with respect to , and denote derivatives with respect to . The over-bars signify that the quantities are defined at the cosmological background level. The trace of Eqn. (2) becomes a dynamical equation for :
(5) |
Recently, a phenomenologically viable model of was proposed in Hu and Sawicki (2007a, b), in which
(6) |
where , and are free dimensionless parameters and
(7) |
This form of can be tuned to match the CDM expansion history while not introducing a true cosmological constant. In addition, the above form has the potential to evade solar system tests via the so-called chameleon mechanism (Khoury and Weltman, 2004a, b), thus remaining viable even with the stringent constraints placed by such tests. However, the HS model has interesting deviations from CDM in the intermediate scales that result in distinct signatures in large scale structure formation. Such deviations have been studied in the linear regime (Song et al., 2007a, b; Pogosian and Silvestri, 2008), but the lack of insights into the modification’s behavior in the non-linear regime places a strong limitation on the applicability of comparison of such studies to observations.
The free parameters, , , and , are chosen so that the resulting background expansion history closely matches to that of CDM. This requirement is equivalent to requiring the background field strength today, , to be sufficiently small () or, alternatively, . In this parameter regime, the field is always near the minimum of the effective potential,
(8) |
In the limit of CDM, we have
(9) |
and thus, to approximate the CDM expansion history, we set
(10) |
Two free parameters remain, namely and , that control how closely the expansion history matches that of CDM. From the high curvature limit of ,
(11) |
we see that larger results in closer match to CDM (i.e., approaches a constant as ) and smaller mimics CDM until later in the expansion history. It is also important to note that the condition requires that at the background level throughout the expansion history because is a monotonically decreasing function of time. Thus, in our numerical simulations, we are able to simplify the field equations given the high curvature assumption.
Phenomenologically viable models require that be small so as to not significantly alter the cosmic expansion history, limiting the field strength to with current observations (Song et al., 2007b) and to with future weak lensing surveys (Hu and Sawicki, 2007a). An even stronger constraint can be placed on the field by considering the dynamics of our galaxy, resulting in an upper bound of (Hu and Sawicki, 2007a). However, such analyses do not self-consistently model large scale clustering that potentially relaxes the tight upper bound.
In the high curvature regime, our model can be expanded as,
(12) | |||||
valid when . Furthermore, the field equation (5) can be simplified by neglecting the small term, thus resulting in
(13) |
Assuming the Friedmann-Robertson-Walker (FRW) metric, we subtract the spatially constant background quantities and thus the field equation becomes an equation for the perturbations,
(14) |
where , , , and the term proportional to has been neglected because . Finally, we work in the quasi-static limit in which , and the field equation is reduced to
(15) |
The field is coupled to the gravitational potential, , via
(16) |
for which the derivation can be found in Hu and Sawicki (2007a). Equations (15) and (16) define a system of equations for the gravitational potential, the solution to which can be used to evolve matter particles. Note that in the limit of ordinary GR, , and the equation for the gravitational potential reduces to the unmodified equation,
(17) |
The strategy that we use to simulate cosmology boils down to solving the two field equations, Eq. 15 and 16, and using the resulting peculiar potential to update the dark matter density.
Iii Numerical Methods
Our N-body simulation is based on the Particle-Mesh algorithm (hereafter PM), in which the particle densities are interpolated onto a regular grid and the potential is solved only at the grid points (see, e.g., Hockney and Eastwood, 1981; Bertschinger, 1998). The choice to use a field representation of the potential instead of particle pair representation (as in Tree codes) is motivated by the fact that the change in the inter-particle forces are mediated by an auxiliary scalar field, namely . The force law between a pair of particles in the presence of depends on the entire behavior of the field between the two particles and thus cannot be expressed as a simple function of their separation. In this section, we outline the method we developed to obtain by numerically solving equations (15) and (16).
iii.1 Code Units
The code units are based on the convention used in Shandarin (1980); Kravtsov et al. (1997), where
and
(18) | |||||
(19) | |||||
(20) | |||||
(21) | |||||
(22) | |||||
(23) |
In the above definitions, is the comoving simulation box size in Mpc, is the number of grid cells in each direction, is the Hubble parameter today, is the critical density today, and is the fraction of non-relativistic matter today relative to the critical density.
Before converting the field equations to code units, we restore the factors of in Eqn (15), which results in
(24) |
Substituting all physical quantities and derivatives in Eqn (24) by their code unit equivalent yields
(25) |
where
(26) |
The quantity is the speed of light in code units, i.e., . Similarly, we transform Eq. (16) into code units and obtain
(27) |
From this point on, we work in code units unless otherwise specified, and we drop the tilde from code variables.
iii.2 Solving for and
The field equation for is a non-linear elliptical partial differential equation because the function is non-linear in general. Because of this non-linearity, standard spectral methods (e.g., FFT method) are not applicable and we therefore resort to a non-linear relaxation method for the solution. Given a general partial differential equation,
(28) |
where is a non-linear differential operator on the field and is the source term, let us define
(29) |
The problem of finding the solution, , for Eq. (28) becomes a problem of finding the root of . The root-finding is carried out by an iterative Newton-Raphson algorithm, which defines the relaxation method (see, e.g., Press et al., 1992),
(30) |
The only remaining difficulty is determining the appropriate form of and its discretization. For our model of , the differential operator is defined as
(31) |
A naive and straightforward discretization of Eq. 25 suffers from severe numerical instability. The differential operator in Eq. 31 is not defined for non-negative values of (see Eq. 12), and thus the iterative solution fails as soon as an iteration oversteps the true solution into the forbidden region. We find that placing an artificial ceiling on the value of does not adequately solve the problem because the stability becomes conditional on very precise tuning of the ceiling due to the steepness of near . The best way to avoid the problem in this model is to redefine the field as
(32) |
and solve for instead of . With this substitution, the new field has infinite domain and does not suffer from abrupt instabilities. The differential operator thus becomes
(33) |
In our implementation, the differential operator is discretized as a variable coefficient Poisson equation, in which the value of at a grid cell identified by , and is updated using the red-black Newton-Gauss-Seidel (NGS) scheme (Press et al., 1992). The actual discretization is too long to be displayed in the body of the text, and is instead placed in Appendix A. We use periodic boundary conditions to close the NGS scheme at the edges of our simulation boxes.
One potential problem with this substitution is the fact that we added additional non-linearity into a problem that is already highly non-linear. In particular, the additional non-linearity is exponential, and thus the errors in the discretized solution, , are potentially enhanced exponentially. In our experience running cosmological simulations with modification, we find that the effects of this addition are negligible in all cases and the iterations converge to the same solution whether this substitution is made or not, while avoiding the singularity of the curvature field at .
Given the solution to using the NGS algorithm and the density field, we can compute the solution to using Eq. 27. The equation is a standard linear Poisson equation and thus is solved efficiently using the fast Fourier transform method. We adopt the cloud-in-cell density interpolation algorithm (hereafter CIC), which treats a particle as a uniform density cube that is equal in size to the underlying grid cell, to compute the density field.
iii.3 Multigrid Acceleration
In general, relaxation methods for solutions of partial differential equations are very inefficient and computationally expensive. However, this inefficiency can be mostly alleviated by the use of multigrid techniques. The multigrid method (Brandt, 1973) is a class of algorithms to solve partial differential equations using a hierarchy of discretizations with increasingly coarse grids. In particular, it exploits the fact that a typical relaxation schemes, such as the NGS scheme previously detailed, smooth out the high frequency residual errors faster than the low frequency error modes. Thus, by interpolating the low frequency error modes onto a coarser grid, the same error mode now appears as a high frequency mode on the coarser grid and can be efficiently reduced by relaxation passes on that grid.
In our working code, we implement a particular type of multigrid method called the Full Approximation Scheme (Briggs et al., 2000, hereafter FAS), which is particularly suited for non-linear boundary value partial differential equations. In order to define the FAS algorithm, we first define two grids, the fine grid (), and the coarse grid (). In addition, we define two inter-grid operators, one that restricts a fine grid field onto a coarse grid, and one that interpolates a coarse grid field onto a fine grid. We label the former as and the latter as , and leave the details of their implementations to Appendix B. Finally, we denote the quantities defined on a coarse grid by a superscript and the quantities defined on a fine grid by a superscript .
The two grid FAS scheme can be summarized by the following. Given the initial guess to the solution, , the FAS algorithm can be summarized as follows:
In a true multigrid method, the third step above is also carried out using a multigrid method, thus forming the recursive hierarchy of coarser grids. The recursion is stopped when the number of coarse grid cells per dimension reaches four; the solution at this level is obtained by iterating NSG until convergence. All relaxations are carried out using the NGS scheme that is appropriate for the given coarse (or fine) grid.
In our implementation, we use the full-weighted restriction operator () and its adjoint, the bilinear interpolation as the prolongation operator (). In practice, we find that the choice of the restriction and prolongation operators does not affect the stability of convergence for most typical cosmological simulation. However, we find the convergence rate when using full-weighting and bilinear interpolation to be better than other, simpler choices.
iii.4 Particle Dynamics
Given the solution for field, the Newtonian potential is computed by solving Eqn (27), which is a linear Poisson equation with a non-standard (i.e., non-GR) source term. Because of its linearity, such equation can be solved efficiently using a fast Fourier transform method (FFT, see Hockney and Eastwood, 1981). The forces on the dark matter particles are computed using Newton’s equations in comoving code coordinates,
(34) | |||||
(35) |
Because the expansion history of our model is tuned to that of flat CDM, is computed simply as
(36) |
where the over-dot denotes a derivative with respect to coordinate time in code units. The gradient of the gravitational potential is computed using a standard finite difference scheme consistent with our grid density assignment operator, namely the CIC scheme. We evolve the particles using a second-order accurate leap-frog method (Klypin and Holtzman, 1997).
Iv Tests of the code
Our N-body simulation is implemented in C++ with shared memory parallelization using OpenMP. In the following subsections, we present various tests of the implementation to assess the correctness of the code as well as its computational speed.
iv.1 Analytically tractable case
First, we construct a density field that has a known solution to check that the code can recover the true solution. In 1D, the field equation (25) becomes
(37) |
which, after setting all constants to 1 and taking and , is
(38) |
We take the analytic solution to be
(39) |
and substitute the solution into equation (38) to find
(40) | |||||
Thus, setting the density field to the expression in (40) and numerically solving for allows us to compare the code accuracy against the known solution, Eq. 39.
In Fig. 1, we show the results of applying our code to this problem with three different grid resolutions. In the left panel, fractional discretization errors, defined as , are shown, and the residual errors, defined as , are shown in the right panel. The average discretization error decreases by approximately an order of magnitude for each two-fold increase in grid resolution. Thus, for , we expect the discretization error to be on the order of , and therefore gives us a useful stopping criteria for solution convergence based on the requirement that the discretization error is the dominant source of error. We terminate the iterative multigrid process when the residual error (which can be computed without knowing the true solution) is smaller than , guaranteeing the residual error to be much smaller than discretization error for our highest resolution simulations.
iv.2 Point mass
In this test, we place a point mass (consistent with the CIC scheme) at the origin of a periodic simulation box, with the corresponding density field given by
(41) |
where the indices , , and , refer to the grid cell indices in , , and directions, respectively. In the subsequent discussion in this section, we restore the tilde notation to distinguish code variables and physical variables. With this density configuration, outside the grid cell at the origin, we expect that
(42) | |||||
where is the effective field mass of the field, given in physical coordinates by
(43) |
for small . The quantity is the second derivative of with respect to . Thus, the solution to Eq. (42) is a Yukawa like potential,
(44) |
where
(45) |
is the effective comoving field mass in code units.
Fig. 2 shows the field solution computed on a grid. The simulation parameters were: Mpc, , , , and . On the same plot, the analytic solution with the expected is plotted. The discrepancy at Mpc is due to the fact that the approximation is only valid for small and thus is only appropriate at large distances away from the central near point mass. Hence, the Yukawa solution in Fig. 2 is normalized such that the solution best matches the simulation at Mpc where the becomes sufficiently small. At Mpc, the solution approaches the level of discretization error for () and hence begins to show deviations from the analytic solution.
iv.3 Dynamics
As previously described in §III.4, the particle trajectories are integrated using a second-order accurate leap-frog scheme. In Fig. 3, we show the orbits of test particles orbiting around a central point mass. The test particles are placed at two and seven grid cells away from the point mass and are given initial velocity such that they would be in circular orbits if there were no integration errors and no modifications. The solid lines show the simulated orbit for the two different radii. Far away from the point source, where the forces are adequately resolved, the orbit is nearly circular as expected for unmodified GR. The smaller orbit, with a theoretical radius of two grid cells, is not well resolved by the underlying potential grid. However, the simulated orbit is still nearly circular, with only a few percent deviation after two complete orbits. The two-body force law, and its modifications due to models, will be further studied in §V.1.
Another popular test of the particles dynamics as well as of the accuracy of the potential solver is the Zeldovich pancake collapse test (Klypin and Shandarin, 1983; Efstathiou et al., 1985; Kravtsov et al., 1997). The evolution of a 1D sine wave is described exactly by the Zeldovich solution,
(46) | |||||
(47) |
where , is the linear growth function, is the time step used in the simulation, and is a normalization constant that fixes the time at shell crossing. By setting the initial conditions at according to the Zeldovich solution, we can compare the subsequent particle evolution to the same solution to gauge the accuracy of our particle dynamics.
For our test, we use , , , and no modifications. The simulation box is set to Mpc, and we use particles and grid cells. Fig. 4 shows the particle momenta versus the particle positions at the time of shell crossing, which is set to be . The dots denote the simulation output and the line represents the exact solution given by the Zeldovich approximation. As expected, the agreement between our code and the Zeldovich solution is good, with maximum deviation less than 1%.
As another check of the particle dynamics, we have run several cosmological simulations without modifications to check our results against the Smith et. al. fits to the power spectrum (Smith et al., 2003). We run three simulations with box sizes () of 256 Mpc (L256), 128 Mpc (L128), and 64 Mpc (L64). All simulations have the following cosmological parameters resembling the third year Wilkinson Microwave Anisotropy Probe (Spergel et al., 2007, WMAP3) results: , , , km/s/Mpc, and . We take the slope of the primordial power spectrum, , to be 0.958.
The initial conditions for the simulations are created using Enzo (O’Shea et al., 2004), a publicly available cosmological N-body + hydrodynamics code. Enzo uses the Zeldovich approximation to displace particles on a uniform grid according to a given initial power spectrum. We use the initial power spectra given by Eistenstein & Hu (Eisenstein and Hu, 1998), not including the effects of baryon acoustic oscillations. The simulations are started at .
All simulations are run with 512 grid cells in each direction (i.e., ) and with particles. Thus, the formal spatial resolutions of the simulations are 0.5 Mpc, 0.25 Mpc, and 0.125 Mpc for the largest, middle, and smallest boxes, respectively. The corresponding mass resolutions are , , and .
For each simulation box configuration, we run multiple simulations with different realization of the initial power spectrum in order to reduce finite sample variance. Two L256 simulations, three L128 simulations, and four L64 simulations were run.
In Fig. 5, we show the power spectrum of our simulations as well as the Smith et. al. fit. Multiple simulations of the same box size are averaged in the figure. The vertical lines denote the half-Nyquist scale of the particles, plotted as a rough guide to show the resolution limits of our simulations. We see from the plot that our simulations match the Smith et. al. fit to within % at most scales smaller than the half-Nyquist scale. The large scatter at low is seen because only a small number of the largest wavelength modes can fit in a simulation box.
iv.4 Time derivative
In our simulation, we work in the quasi-static limit of the full field equation for , Eq. 5, in which we neglect the time derivative term. There is no simple way to check the validity of our assumption short of running a full simulation including the time derivative term; nevertheless, we can check that our assumption is at least self-consistent by running modified cosmological simulations with quasi-static assumption and comparing the time derivative of the field to its spatial derivative.
We use the same simulation setup as the cosmological simulations in §IV.3 except for the inclusion of modification with . The time derivative of the field is estimated by finite differencing of fields at successive time steps. In Fig. 6, we show the fraction of the 2-norms of the time derivative and spatial derivatives of the field at different scale factors. Although the ratio varies by time, we see that the time derivative is generally much smaller than than the spatial derivative and thus that our simulations are consistent with our assumptions.
iv.5 Code timing and memory requirement
Even with the efficiency of the multigrid method, solving a 3D non-linear partial differential equation with reasonable resolution is still a computationally expensive endeavor. In Fig. 7, we plot the residual error as a function of the CPU time spent during the NGS and multigrid iterations of the point mass solution studied in §IV.2. The computations were carried out on a 2.8 GHz AMD Opteron 2220 processor without using OpenMP multi-threading (i.e., single threaded). For the NGS method, the high frequency error modes are reduced more efficiently than low frequency modes. Thus, the residual error goes down fast during the initial few hundred seconds of the NGS iterations, until the low frequency error modes dominate the overall residual error budget. Thereafter, the iterations converge slowly, rendering the method unacceptable for cosmological simulations in which the field must be evaluated hundreds of times. With the multigrid method, because of its use of coarse grids, error modes of all frequencies decay at nearly the same rate, and thus the overall convergence is greatly accelerated. Improved convergence rates using multigrid method, sometimes by as much as two orders of magnitudes, makes cosmological simulations a reality.
As shown in Fig. 7, convergence to a reasonable error tolerance takes s on an AMD Opteron 2220 processor even for a modest resolution. The FAS scheme has one of the best computation time scaling, growing linearly to the number of grid points (FFT, for comparison, scales as ), and thus a simulation will converge in 10,000 s. The field equation must be solved for each individual time step of the simulation, and thus resulting in potentially thousands of multigrid solutions totaling a few million seconds of CPU time. With our discretization scheme, we find that the simulation code spends of CPU time solving for the field. Such percentage is at best a rough estimate, as the number of multigrid iterations required for convergence depends on the setup of the simulations.
Fortunately, several methods to speed up the simulation are possible. First, the particle positions between successive time steps does not change greatly, and thus the solution should not change significantly as well. The field solution from the previous time step provides a good initial guess for the iterative solver and we find that it results in faster convergence. Secondly, the FAS scheme is easily parallelized using shared memory parallelization with almost perfect scaling with the number of processor cores available. The NGS smoothing used extensively within the FAS scheme can be arranged in the so called red-black ordering (Briggs et al., 2000), such that smoothing of one particular grid point is independent of all other points except for its six nearest neighbors. This locality allows the smoothing to be split into local domains, each of which can be computed on separate processors. Combining such acceleration methods, a typical simulation can be carried out in roughly 10 days on our hardware. We see no difficulties in extending the parallelism to distributed memory, although we have not implemented such extension.
Memory requirement for the FAS scheme can be high compared to traditional PM simulations. The hierarchy of coarser grids add to the storage requirement for the field compared to the normal storage requirement for a 3D grid scaling as . Furthermore, in order to accelerate the computation, we store the current residual error and use one additional grid hierarchy of the same size for transient computations. For a grid, the required memory for our FAS implementation is GiB using double precision floating point format (64 bits each). For comparison, storage of the density field requires 1 GiB, and the storage of dark matter particles requires 6 GiB for particles.
V Results
In this section, we study a few simulations in order to highlight the differences between and ordinary CDM physics. A more comprehensive study of the non-linear physics is detailed in a companion paper (Oyaizu et al., 2008).
v.1 Two-body force laws
The first interesting case is the study of two body force laws in cosmology. As is well known, in ordinary GR, two-body forces decay as as the separation between the bodies, , is varied. However, in cosmology, the new field mediates a fifth force that modifies two-body interactions and potentially change the large scale structures of the Universe.
Such deviations from GR are easy to see using the same setup as the near point mass text in § IV.2. We solve for the peculiar gravitational potential, , with the field and compute the gravitational acceleration at random positions in the simulation box, labeling them . In addition, we compute the standard Newtonian force using the inverse square law for the same set of points, labeling them . We use a 400 Mpc computation box with , , and . Fig. 8 shows the relative differences between and as a function of the distance away from the point mass. For strong fields, the gravitational forces near the central point mass are enhanced by as much as over the unmodified forces. The range of the force enhancement is determined by the Compton wavelength of the field, which depends on the local field mass as
(48) |
Thus, the stronger the field, the farther out the gravitational enhancement of the field propagates, which is confirmed by Fig. 8. For the very weak field with , the Compton wavelength is much smaller than the computational grid cells and therefore the gravitationally enhanced regions are not resolved.
In the strong field limit, the large Compton wavelength of the field makes it smooth and stiff, resisting the fluctuations in the underlying density field. Thus, in this limit, the deviation of the field, and hence of the curvature field , vanishes and we have . The equation for the gravitational potential, Eq. 16, reduces to
(49) |
We see that the result is an effective enhancement of the Newton’s constant by 33%, and hence the two-body forces enhanced by the same amount. Our two-body numerical solution captures this effect well, with the acceleration curve plateauing at the correct value.
Features not due to the presence of the field are also seen in Fig. 8. The large scatter near are due to the lack of force resolution at scales comparable to the underlying grid cells. At scales even smaller than the grid cells, computed forces vanish; this force decay is essential for the N-body system to remain collision-free. On the other end of the scale, the periodic boundary condition begins to affect the forces as the test masses becomes gravitationally attracted to the neighboring point mass.
At first, Fig. 8 seems to contradict the stringent solar system tests carried out so far and would constrain to very small values. However, the particular model chosen in this study shows the so called chameleon behavior (Khoury and Weltman, 2004a, b), in which the field strength is pushed to near zero in high density regions. Small field strength in the high density region causes the effective field mass to increase and therefore the Compton wavelength to decrease. Thus, in our solar system, which sits in a relatively high density region of the Universe, the chameleon behavior may cause the field to be undetectable using our current technologies.
Our code is expected to reproduce this non-linear chameleon behavior that is not incorporated into the linear analyses carried out so far. In order to investigate such effects, we place a spherical top-hat of constant density at the center of a 200 Mpc simulation box. The sphere is 20 grid cells in radius, corresponding to 15.6 Mpc, and has an uniform over-density of . Random test masses are placed both inside and outside the sphere, and the accelerations of the masses are calculated using the exact solution () and simulations with various field strengths (). The exact solution of the acceleration increases linearly within the sphere and decreases as outside the sphere.
In Fig. 9, we show the relative difference between and as a function of the separation between the center of the spherical blob and the test mass. Outside of the blob, the gravitational force is enhanced as expected from the results of Fig. 8. However, inside the blob, the field quickly dies down and the gravitational forces converge to that of CDM. Strong fields can penetrate the density field farther, as seen in the case of that decays in Mpc.
v.2 Revisiting the Zeldovich pancake
The Zeldovich 1D plane wave collapse offers insights into the differences in the particle dynamics in the presence of modifications. We run three Zeldovich pancake simulations starting with the same cosmology and initial conditions as in Fig. 4 except for modifications. In Fig. 10, we show the resulting difference in particle momenta relative to those of unmodified GR simulation. Because of the fifth force due to , the sine wave collapses faster as the field strength is increased. In addition, we see the effects of the finite range of the fifth force, which is determined by the Compton wavelength of the field, Eq. 48. In the case of weak field, , there is no significant enhancement in the particle momenta at distances greater than code units () away from the collapsed pancake, and the range of the fifth force is shown to be larger for stronger fields. For the strongest field, the fifth force does not decay sufficiently before the periodic image of the pancake starts to become important.
The Zeldovich collapses show no sign of the chameleon behavior, even for the weakest field. The pancake is essentially a sheet mass and thus does not have an extended region of high density. Therefore, the field does not appreciably decay near the pancake and no chameleon behavior is expected.
v.3 Cosmological Power Spectrum
Finally, we apply our code to a realistic cosmological simulation in order to study the effects of the field to the large scale structures. We have run cosmological simulations using the same setup as those found in §IV.3 and with . The strength of the field is chosen to highlight the differences in the power spectra.
With our simulation setup having twice as many grid points as particles per one dimension, we find that at the initial redshift, some of the grid cells in the simulations have no particles nearby and therefore have zero density. At early redshifts, the cosmology is very nearly CDM, and we have . Thus, in the empty grid cells, the curvature is almost exactly zero and we have a formal divergence of the field in such grid cells. The particle dynamics, however, is governed by the peculiar potential, and therefore by the curvature and density fields, not directly by the field (see Eq. 27). We have checked that our potential solver and particle trajectory integrators remain stable and accurate even when such divergences are present (see Fig. 3). At late times, we also find empty grid cells within large voids. However, the mean Compton wavelength of the field becomes sufficiently large and smoothes out the field, thereby avoiding any potential divergences in the field at late times.
Fig. 11 shows the relative differences between the power spectra with and without modifications, along with the linear prediction for the power spectrum enhancement taken from Hu and Sawicki (2007a). The linear prediction assumes no chameleon behavior and thus is expected to overestimate the power spectrum enhancement when there is strong suppression of field in high density regions. The simulated power spectra match the linear predictions at large scales, where non-linear effects, such as the chameleon mechanism, is expected to be negligible. The smallest simulation box, Mpc, does not have sufficient statistics to accurately represent the largest modes. At scales smaller than Mpc, we see a striking deviation of the simulations from the linear prediction, showing clear signs the previously unquantified effects of non-linear physics.
More detailed comparison of the power spectra are carried out in the companion paper.
Vi conclusion
In this paper, we have introduced a method to construct a cosmological N-body simulation with a class of viable modifications to general relativity. Our simulation code, implemented in standard C++, solves the HS model of modification using particle mesh method, particularly suited because the scalar degree of freedom in the model is well represented as a scalar field (i.e., a mesh) coupled to the Newtonian potential. The non-linear field equation for is solved on the mesh using an efficient multigrid scheme that is shown to converge robustly to the desired solution. The code is also shown to reproduce the chameleon behavior of the HS model.
The code should allow us to gain insight into the quasi-nonliear physics of modified gravity models. In general, gravitational forces are enhanced under the HS model, and we show that the resulting power spectra are also enhanced compared to the CDM power spectrum. We further show hints that the power spectrum enhancement is smaller than the linear predictions due to the chameleon mechanism that has not so far been modelled. In a companion paper (Oyaizu et al., 2008), we further study the effects of the HS model on the dark matter power spectrum and its cosmological implications. In addition, effects on the cluster mass function and weak lensing statistics are potentially interesting and should be investigated in the future.
Lastly, the work presented here can be used as a general framework for non-linear simulation of any minimally coupled scalar field model of gravity, provided that the resulting field equations are sufficiently stable under relaxation and multigrid schemes. In particular, an interesting exercise might be to apply this framework to the Dvali, Gabadadze, and Porrati (DGP) model (Dvali et al., 2000).
Vii Acknowledgments
We would like to thank Wayne Hu, Marcos Lima, Justin Khoury, Josh Frieman, and Jeremy Tinker for insightful and useful discussions, and Andrey Kravtsov for guiding the authors during the implementation and testing of the simulation code. This work was supported in part by the Kavli Institute for Cosmological Physics (KICP) at the University of Chicago through grants NSF PHY-0114422 and NSF PHY-0551142 and an endowment from the Kavli Foundation and its founder Fred Kavli. HO was additionally supported by the NSF grants AST-0239759, AST-0507666, and AST-0708154 at the University of Chicago. Some of the computations used in this work have been performed on the Joint Fermilab - KICP Supercomputing Cluster, supported by grants from Fermilab, Kavli Insititute for Cosmological Physics, and the University of Chicago. This research has made use of NASA’s Astrophysics Data System and arXiv.org preprint service funded by Cornell University and the NSF.
Appendix A Details of the discretization
In our implementation, Eq. (33) is discretized much like a standard variable coefficient Poisson equation. We let , , and . Then, we have
(50) | |||||
We find that the above discretization performs better than discretizing the equation after expanding out the derivative operators. Even using this discretization, we find that the iterations can diverge if the initial guess is sufficiently poor. We find that using variable, and often smaller, Newton-Raphson step sizes (i.e., successive under-relaxation) can effectively avoid such divergences until the trial solution is sufficiently close to the true solution.
Appendix B FAS implementation
In our FAS implementation, we use the full-weighting restriction operator () and the bilinear interpolation operator () for prolongation. The full-weighting operator, in 3D, is given by
(51) |
where
Note that the superscript denotes a variable defined on the coarse grid (i.e., a grid with spacing of ) and the superscript denotes a variable defined on the fine grid. The bilinear interpolation operator is given by
Other choices of restriction and prolongation operators are certainly possible. However, we find that the above choices result in the fastest and most stable FAS implementation for this particular problem.
References
- Frieman et al. (2008) J. Frieman, M. Turner, and D. Huterer, ArXiv e-prints 803 (2008), eprint 0803.0982.
- Carroll et al. (2004) S. M. Carroll, V. Duvvuri, M. Trodden, and M. S. Turner, Phys. Rev. D 70, 043528 (2004), eprint arXiv:astro-ph/0306438.
- Nojiri and Odintsov (2003) S. Nojiri and S. D. Odintsov, Phys. Rev. D 68, 123512 (2003), eprint arXiv:hep-th/0307288.
- Capozziello et al. (2003) S. Capozziello, S. Carloni, and A. Troisi, ArXiv Astrophysics e-prints (2003), eprint astro-ph/0303041.
- Williams et al. (1996) J. G. Williams, X. X. Newhall, and J. O. Dickey, Phys. Rev. D 53, 6730 (1996).
- Anderson et al. (1998) J. D. Anderson, P. A. Laing, E. L. Lau, A. S. Liu, M. M. Nieto, and S. G. Turyshev, Physical Review Letters 81, 2858 (1998), eprint arXiv:gr-qc/9808081.
- Bertotti et al. (2003) B. Bertotti, L. Iess, and P. Tortora, Nature (London) 425, 374 (2003).
- Williams et al. (2004) J. G. Williams, S. G. Turyshev, and D. H. Boggs, Physical Review Letters 93, 261101 (2004), eprint arXiv:gr-qc/0411113.
- Knop et al. (2003) R. A. Knop, G. Aldering, R. Amanullah, P. Astier, G. Blanc, M. S. Burns, A. Conley, S. E. Deustua, M. Doi, R. Ellis, et al., Astrophys. J. 598, 102 (2003), eprint arXiv:astro-ph/0309368.
- Riess et al. (2004) A. G. Riess, L.-G. Strolger, J. Tonry, S. Casertano, H. C. Ferguson, B. Mobasher, P. Challis, A. V. Filippenko, S. Jha, W. Li, et al., Astrophys. J. 607, 665 (2004), eprint arXiv:astro-ph/0402512.
- Spergel et al. (2007) D. N. Spergel, R. Bean, O. Doré, M. R. Nolta, C. L. Bennett, J. Dunkley, G. Hinshaw, N. Jarosik, E. Komatsu, L. Page, et al., Astrophys. J. Supplement 170, 377 (2007), eprint arXiv:astro-ph/0603449.
- Hu and Sawicki (2007a) W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007a), eprint arXiv:0705.1158.
- Khoury and Weltman (2004a) J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 (2004a), eprint arXiv:astro-ph/0309411.
- Khoury and Weltman (2004b) J. Khoury and A. Weltman, Physical Review Letters 93, 171104 (2004b), eprint arXiv:astro-ph/0309300.
- Song et al. (2007a) Y.-S. Song, W. Hu, and I. Sawicki, Phys. Rev. D 75, 044004 (2007a), eprint arXiv:astro-ph/0610532.
- Pogosian and Silvestri (2008) L. Pogosian and A. Silvestri, Phys. Rev. D 77, 023503 (2008), eprint arXiv:0709.0296.
- Hu and Sawicki (2007b) W. Hu and I. Sawicki, Phys. Rev. D 76, 104043 (2007b), eprint arXiv:0708.1190.
- Song et al. (2007b) Y.-S. Song, H. Peiris, and W. Hu, Phys. Rev. D 76, 063517 (2007b), eprint arXiv:0706.2399.
- Hockney and Eastwood (1981) R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (Computer Simulation Using Particles, New York: McGraw-Hill, 1981, 1981).
- Bertschinger (1998) E. Bertschinger, Annual Review of Astronomy and Astrophysics 36, 599 (1998).
- Shandarin (1980) S. F. Shandarin, Astrophysics 16, 439 (1980).
- Kravtsov et al. (1997) A. V. Kravtsov, A. A. Klypin, and A. M. Khokhlov, Astrophys. J. Supplement 111, 73 (1997), eprint arXiv:astro-ph/9701195.
- Press et al. (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes in C. The art of scientific computing (Cambridge: University Press, —c1992, 2nd ed., 1992).
- Brandt (1973) A. Brandt, Proceedings of Third International Conference on Numerical Methods in Fluid Mechanics (1973).
- Briggs et al. (2000) W. L. Briggs, V. E. Henson, and S. F. McCormick, A multigrid tutorial (2nd ed.) (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2000), ISBN 0-89871-462-1.
- Klypin and Holtzman (1997) A. Klypin and J. Holtzman, ArXiv Astrophysics e-prints (1997), eprint astro-ph/9712217.
- Klypin and Shandarin (1983) A. A. Klypin and S. F. Shandarin, Mon. Not. R. Astron. Soc. 204, 891 (1983).
- Efstathiou et al. (1985) G. Efstathiou, M. Davis, S. D. M. White, and C. S. Frenk, Astrophys. J. Supplement 57, 241 (1985).
- Smith et al. (2003) R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, Mon. Not. R. Astron. Soc. 341, 1311 (2003), eprint arXiv:astro-ph/0207664.
- O’Shea et al. (2004) B. W. O’Shea, G. Bryan, J. Bordner, M. L. Norman, T. Abel, R. Harkness, and A. Kritsuk, ArXiv Astrophysics e-prints (2004), eprint astro-ph/0403044.
- Eisenstein and Hu (1998) D. J. Eisenstein and W. Hu, Astrophys. J. 496, 605 (1998), eprint arXiv:astro-ph/9709112.
- Oyaizu et al. (2008) H. Oyaizu, M. Lima, and W. Hu, ArXiv e-prints (2008), eprint astro-ph/0807.2462.
- Dvali et al. (2000) G. Dvali, G. Gabadadze, and M. Porrati, Physics Letters B 485, 208 (2000), eprint arXiv:hep-th/0005016.