Einstein-Vlasov Calculations of Structure Formation

TL;DR

This paper compares fully relativistic Einstein-Vlasov simulations with Newtonian N-body simulations to assess the accuracy of Newtonian methods in modeling cosmic structure formation, finding good agreement in most regimes.

Contribution

It provides the first detailed comparison between Einstein-Vlasov and Newtonian N-body simulations for structure formation, highlighting relativistic effects at large inhomogeneities.

Findings

Newtonian and relativistic results agree well at small inhomogeneities.

Large inhomogeneities cause faster collapse in Newtonian simulations.

Outside high potential regions, spacetime evolution differences are minimal.

Abstract

We study the dynamics of small inhomogeneities in an expanding universe collapsing to form bound structures using full solutions of the Einstein-Vlasov (N-body) equations. We compare these to standard Newtonian N-body solutions using quantities defined with respect to fiducial observers in order to bound relativistic effects. We focus on simplified initial conditions containing a limited range of length scales, but vary the inhomogeneities from small magnitude, where the Newtonian and general-relativistic calculations agree quite well, to large magnitude, where the background metric receives an order one correction. For large inhomogeneities, we find that the collapse of overdensities tends to happen faster in Newtonian calculations relative to fully general-relativistic ones. Even in this extreme regime, the differences in the spacetime evolution outside the regions of large gravitational potential and velocity are small. For standard cosmological values, we corroborate the robustness of Newtonian N-body simulations to model large scale perturbations and the related cosmic variance in the local expansion rate.

Equations

$${\delta }_N=\sum _i{\bar{\delta }}_i\sin (k{x}^i),$$

$$v^i=H_0{\delta }_i\cos (k{x}^i)/k.$$

$$\delta x^i=-\frac{4\pi }{{\rho }_0}{\partial }_i{\Psi }_N(a=1),$$

$${\partial }^i{\partial }_i{\Psi }_N=4\pi a^2{\rho }_0{\delta }_N,$$

$${\delta }_{\mathrm{obs}}(\tau ):=(\rho (\tau )/{\rho }_0)a_p^{-3}-1,$$

$$a_p:={\left[3\tau H_0/2+1\right]}^{2/3}$$

$$z=-1+\frac{(u_a{k}^a{)}_{\mathrm{emit}}}{(u_a{k}^a{)}_{\mathrm{obs}}}.$$

Full text

I Introduction

Recently, there has been a growing interest in quantifying the importance of effects that are both nonlinear and relativistic on the large scale evolution and development of structure in the Universe Bentivegna and Bruni (2016); Giblin et al. (2016a); Rekier et al. (2015); Giblin et al. (2016b); Macpherson et al. (2017); Daverio et al. (2016); East et al. (2018); Macpherson et al. (2019); Giblin et al. (2019); Daverio et al. (2019). This means studying effects that may be missed by the standard tool for studying cosmological structure formation: Newtonian $N$-body simulations. The motivation for such studies ranges from answering claims that small scale nonlinearities may have a strong “backreaction” on large scales on the one extreme Buchert (2000); Kolb et al. (2005); Rasanen (2011); Ishibashi and Wald (2006); Green and Wald (2014), to the desire to quantify small, subpercent relativistic effects which may soon become observable in the era of precision cosmology Laureijs et al. (2011); LSST Dark Energy Science Collaboration (2012); Collaboration (2017).

There are a number of challenges in performing a full, nonperturbative general-relativistic (GR) calculation of structure formation. Solving the Einstein equations requires both solving a set of constraint equations (typically elliptic) at the initial time and evolving hyperbolic equations for the metric which have characteristics that propagate at the speed of light. The latter imposes a severe restriction on the timestep of the simulation compared to the case where the gravity is completely determined by an elliptic equation and the matter moves nonrelativistically. Resolving the small scales of collapsed structures is already very challenging within the Newtonian framework Heitmann et al. (2008); Schneider et al. (2016), and this restriction makes the GR case much more severe. Hence, most calculations beginning with a range of length scales very quickly become underresolved. One approach is to only include some general-relativistic corrections which do not break the elliptic description of gravity Adamek et al. (2016); Barrera-Hinojosa and Li (2019). However, this requires making a priori assumptions about which terms can be neglected.

GR simulations also tend to discretize the metric functions on grids, which makes it natural to use a fluid description of the cold dark matter which can be discretized on the same grid. This is what has been done for most full GR calculations of cosmological structure to date (Refs. Pretorius and East (2018); Giblin et al. (2019); Daverio et al. (2019) are exceptions to this). However, such fluid descriptions break down as soon as multistream regions emerge, which of course are generic features of structure formation.

Finally, there is the difficulty of distinguishing and quantifying the magnitude of effects coming from nonlinear gravity, from those solely due to nonlinear perturbations in the matter (which will be captured by standard Newtonian calculations) East et al. (2018). For example, one cannot simply look at how inhomogeneous various functions of the metric are in a GR simulation. Related to this, when one is considering nonlinear deviations from a homogeneous spacetime, coordinate ambiguities make it difficult to interpret the metric functions directly, and one has to be careful to compute gauge invariant quantities in order to make a meaningful comparison Giblin et al. (2019).

This work extends that of Ref. East et al. (2018), where a direct comparison of Newtonian and GR simulations of structure formation was performed utilizing the dictionary of Refs. Chisari and Zaldarriaga (2011); Green and Wald (2012) to generate consistent initial conditions in both simulations and to compare observables. In Ref. East et al. (2018), a fluid description of the matter was used for the GR calculations, which meant that the comparison became unreliable past the point where multistream regions would develop. Here, we use the methods of Ref. Pretorius and East (2018) to solve the Einstein-Vlasov equations, allowing us to continue the comparison as bound structures are formed. We sidestep some of the computational challenges mentioned above by considering simplified initial conditions, where the perturbations are concentrated at a single wavelength, but consider various magnitudes for the inhomogeneities. For large enough inhomogeneities (in excess of standard cosmological values), we do find appreciable deviations between the Newtonian and GR calculations, with the collapse of overdensities happening faster in the former. However, in the regime where this occurs, it is already clear from the Newtonian calculation itself that deviations are expected since the gravitational potential and velocities relative to the speed of light are becoming comparable to unity. Furthermore, even in such cases, we find that outside the regions of large gravitational potential, the agreement between the two methods in observables like the evolution of the density and the propagation of light is still good.

The remainder of this paper is as follows. In Sec. II, we describe the initial conditions we consider, the methods we use to evolve in both a full GR and Newtonian framework, and the diagnostic quantities we use to compare the two. In Sec. III, we present the results of our calculations evolving inhomogeneities of various magnitudes, and in Sec. IV we conclude. In the appendix, we present results estimating the numerical errors in our calculations. We use units with $G=c=1$ throughout.

II Methodology

II.1 Initial conditions

Following Refs. Bentivegna and Bruni (2016); East et al. (2018), we consider a simple set of initial conditions consisting of density perturbations about a homogeneous solution. The homogeneous solution is characterized by its initial expansion rate $H_{0}$, and hence and density $\rho_{0}:=8\pi/3H_{0}^{2}$, which sets the overall scale. The perturbations are taken to be in each of the Cartesian directions with initial wavelength that is four times the Hubble radius at the beginning of the calculation. That is, we take the Newtonian density contrast to be

[TABLE]

with $k=\pi H_{0}/2$. We introduce a small asymmetry between the different Cartesian directions by letting $\delta_{i}=\bar{\delta}(1,\ 0.9,\ 1.1)$, and we consider varying magnitude density perturbations $\bar{\delta}\times 10^{2}=0.25$, 0.5, 1, and 5. The initial velocity is given by the Zel’dovich approximation (Zel’dovich, 1970)

[TABLE]

These initial conditions have a maximum overdensity at $\left(0,0,0\right)$ and maximum underdensity at $\left(\pi/k,\pi/k,\pi/k\right)$.

As described in detail in Ref. East et al. (2018), fully general-relativistic initial data are calculated using the dictionary of Refs. Chisari and Zaldarriaga (2011); Green and Wald (2012) to determine the approximate metric and stress-energy tensor, and then solving the full Einstein constraint equations in the conformal thin-sandwich formulation East et al. (2012) for any nonlinear corrections.

II.2 Newtonian simulations

The Newtonian $N$-body simulations are performed using the GADGET-2 code Springel (2005) with a TreePM algorithm for the gravity solver (Xu, 1995). These simulations serve as a reference to standard computational cosmology, where the evolution of the cosmic density field is governed by Newtonian gravity, and is fully separated from the background expansion, described in turn by the Friedmann equation. GADGET-2 has been validated in a number of comparison studies verifying the accuracy and robustness of various numerical implementations of cold dark matter cosmological simulations (see, e.g., Refs. Heitmann et al. (2008); Kim et al. (2014); Heitmann et al. (2005)).

We generate conditions by displacing particles from a regular grid according to the field given by the Zel’dovich approximation (Zel’dovich, 1970)

[TABLE]

where $\Psi_{N}$ is the Newtonian gravitational potential given by

[TABLE]

and by convention the scale factor $a$ is set to unity at the beginning of the calculation. The resulting density field that is inferred from the positions of the particles reproduces the input density up to the second order corrections in the density contrast. As in Ref. (East et al., 2018), we apply the corrections by means of a minimal adjustment of particle’s masses. The particle masses are set in such a way that they compensate all local differences between the actual (as calculated by the employed density estimator, described below) and input density evaluated at the position of every particle. We note that the introduced corrections are small (subpercent level), but they guarantee a high-accuracy match between initial conditions of the Newtonian and GR simulations.

The density field is not explicitly evolved in the $N$-body simulations, and it can only be derived from the positions of the particles. Here, we employ a well-tested method for measuring matter density in cosmological simulations of cold dark matter, based on tracing the evolution of the Lagrangian tessellation of the dark matter manifold in phase space (Shandarin et al., 2012; Abel et al., 2012). Density is estimated by means of scaling the initial density according to a relative change of the volume of tetrahedral mass elements defined in the initial tessellation. In single-stream regions (no shell crossing), local density at a given position is determined solely by a single tetrahedral cell containing this point, while density in multistream regions (after shell crossing) arises from multiple density contributions coming from all locally overlapping tetrahedral cells.

The employed density estimator outperforms more traditional techniques such as cloud-in-cell (CIC) in several respects. Here, we emphasize that the estimator can be applied locally, and it does not suffer from undersampling in single-stream regions, making it an ideal method for tracing the density field in voids. On the other hand, density estimates in multiple-stream regions should be regarded with reservation, because the full robustness of the estimator requires simulations with a computationally heavy adaptive refinement of tessellation cells (Hahn and Angulo, 2016). In particular, density estimation in the center of dark mater haloes depends on resolution, and there is no guarantee that the computation can converge due to the cuspy nature of dark matter density profiles, although precision estimated from comparing results based on different tessellations at fixed resolution is of the order of 0.1 dex (Abel et al., 2012). The problem of resolution dependence can be circumvented by employing a density estimator with a fixed smoothing scale in comoving coordinates instead. Bearing this in mind, we include CIC estimates of density in some cases for comparison with the GR calculation (which does not utilize tetrahedral cells).

Unless otherwise stated, the results shown here are obtained used $N=196^{3}$ particles. We also run select cases using $N=128^{3}$ in order to estimate numerical errors. The simulations were carried out with a force softening of $5\times 10^{-4}$ (high resolution) and $8\times 10^{-4}$ (low resolution) in units of the simulation domain length $L$.

In order to compute the trajectories of freely falling test particles, we follow the evolution of the tetrahedral cells containing the initial positions of the test particles. The positions of the evolved test particles are then computed by interpolating between the displacements of cells vertices, which are always given by dark matter particles.

II.3 GR simulations

The fully general-relativistic $N$-body simulations are performed using the methods described in Ref. Pretorius and East (2018). This code was also recently used to follow black hole formation from collisionless matter East (2019). As in the Newtonian simulations, we determine the initial particle positions by starting from a uniform lattice of particles and then displacing each particle slightly according to the Zel’dovich approximation (given by Eq. (31) in Ref. Chisari and Zaldarriaga (2011)). However, there will be a small nonlinear correction to the density field which we will need to apply to the particle distribution. To do this, we use slightly nonuniform masses for the particles, given by rescaling the masses in proportion to the ratio of the desired density to that obtained from the Zel’dovich approximation.

Though the code used here does implement adaptive mesh refinement (see Ref. Pretorius and East (2018)), for this study we restrict to uniform grids. We do this mainly for efficiency, though we note that the results in the appendix indicate that, at late times in our simulations, the numerical error is mainly dominated by the number of particles. For most of the results presented here, we use resolution with 96 points across the wavelength of the initial perturbation, and $4^{3}$ particles per grid cell. However, we run select cases at multiple resolutions utilizing $2/3\times$ and $4/3\times$ as many grid points in order to establish convergence and estimate truncation error. See the appendix for details.

For comparison, we also include a few results that are calculated by treating the matter as a pressureless fluid as described in Ref. East et al. (2018).

II.4 Comparing observables

In order to compare the results of the Newtonian and GR $N$-body evolutions, we compute several quantities defined with respect to fiducial observers, as detailed in  East et al. (2018). We compute the matter density along the worldlines of timelike observers and use this quantity as a function of proper time $\rho(\tau)$ to define an effective density contrast:

[TABLE]

where

[TABLE]

is a convenient parametrization of the proper time using the Lemaître-Friedmann-Robertson-Walker (LFRW) expression for the scale factor that would hold in the homogeneous case. We emphasize that since $H_{0}$ (and hence $\rho_{0}:=8\pi/3H_{0}^{2}$) is just a constant that sets the overall scale of our initial conditions, $\delta_{\rm obs}(a_{p})$ is just a convenient reparameterization of density as a function of proper time.

We also measure properties of the spacetimes using null geodesics which are “emitted” and subsequently “observed” by fiducial timelike observers. If $k^{a}$ is the four momentum of the null geodesic and $u^{a}$ is the four velocity of emitter/observer, we can compute a redshift factor

[TABLE]

For $u^{a}$, we take the four velocity implied by the stress-energy tensor $T^{ab}=\rho u^{a}u^{b}$, which weights the contributions from different particles in the case of multistream regions. We can also use the deviation of neighboring null geodesics to compute the luminosity distance (or, equivalently through the reciprocity relation, the angular distance Etherington (2007)) as a function of the redshift $D_{L}(z)$ along each null ray.

For the GR simulations, these quantities are computed by including extra tracer particles which are evolved in the same way as the matter particles (but without backreacting). For the Newtonian simulations, these quantities are computed by reconstructing the effective spacetime using the Newtonian-GR dictionary of Chisari and Zaldarriaga (2011); Green and Wald (2012) and integrating the resulting geodesic equation. Hence, the Newtonian calculation also includes relativistic effects in the propagation of light, etc., and the comparison is really of how much the spacetimes implied by the two methods of calculation differ.

III Results

With the initial conditions we have chosen, as the spacetime expands and the inhomogeneities move inside the horizon, a growing void emerges at the point of maximum underdensity, and a bound, multistream region (i.e. a halo) is formed at the point of maximum overdensity. In the top and middle panels of Fig. 1, we show the density contrast measure $\delta_{\rm obs}$ at these two points for cases with different magnitudes of the initial inhomogeneities. The Newtonian and GR calculations show good agreement at the underdensity for all cases, even as the density contrast becomes highly nonlinear.

For the overdensity, two differences are noticeable. The first is that the collapse and halo formation occurs slightly earlier for the Newtonian case, and this difference increases as the initial inhomogeneities become larger (and hence more relativistic). The second is that the saturation density is significantly larger for the Newtonian case. We shall not focus too much on the latter since this is fairly sensitive to numerical effects such as the finite number of particles and the smoothing length. In the bottom panel of Fig. 1, we show for the $\bar{\delta}=0.01$ case a comparison of how this quantity changes, both with numerical resolution, and with a particle versus pressureless fluid treatment of the matter. Here it can also be seen that with a CIC estimate of the density, the maximum density contrast for the Newtonian calculation is much closer to the GR result (which similarly deposits each particle’s stress energy on neighboring grid points). In the GR pressureless fluid treatment, the calculation breaks down at shell crossing, whereas with the particle treatment the density eventually saturates. In either case, finite resolution tends to lead an underestimate of the density around this point. However, even taking this into account, the collapse happens faster in the Newtonian case compared to the GR case. This discrepancy increases with increasing inhomogeneity amplitude and becomes quite pronounced for the case with $\bar{\delta}=0.05$. For this extreme case, the Newtonian calculation has to be terminated when the magnitude of the Newtonian potential $\psi$ becomes $\sim 1/2$. We discuss this case in more detail below.

The differences in the evolution of multistream regions can be tracked by considering a set of fiducial observers, comoving with the matter, that are initially displaced from the halo by some distance, and comparing the proper time it takes for them to eventually fall through the point of maximum overdensity and begin to oscillate around it. This is illustrated in Fig. 2, where it is apparent that as the size of the inhomogeneities increases, and the collapse takes place more quickly and at scales more comparable to the Hubble scale, the relative discrepancy between the Newtonian and GR cases increases, with the Newtonian case exhibiting faster collapse. (We note that in general these coordinate distances are gauge dependent, but the time the particles cross the overdensity is not.)

Figure 3 shows the differences between the Newtonian and GR positions of freely falling particles from Fig. 2 as a function of the absolute magnitude of the infall velocity inferred from the Newtonian simulation. For the sake of clarity, we only show the trajectories up until the time where they first cross the halo center in the Newtonian run. The comparison demonstrates that the Newtonian trajectories closely follow their GR counterparts, as long as infall velocities do not exceed the limits of nonrelativistic dynamics. Noticeable discrepancies between the two simulations occur when the particles reach relativistic velocities. The apparent differences reflect the limited accuracy of the Newtonian simulations when there is a violation of the nonrelativistic assumption. Particles in the Newtonian simulations are accelerated to larger velocities, giving rise to a faster collapse onto the central object than in the GR simulations.

We can also compare the differences in the effective spacetimes using the propagation of light. In Fig. 4, we compare the luminosity–redshift relation for fiducial light rays propagating between the points of minimum and maximum density. From the comparison with the homogeneous solution shown in the left column of Fig. 4, one can see that the cases considered here have large, nonlinear deviations from the LFRW behavior. Nevertheless, as evident in the right column, the differences between the GR and Newtonian case remain much smaller, in most cases subpercent and consistent with numerical truncation error (see the appendix and Ref. East et al. (2018)), indicating the differences in the spacetimes are small.

For the larger amplitude inhomogeneities, light rays emitted from the overdensity at later times have a $D_{L}(z)$ that is slightly smaller for the GR calculation than the Newtonian counterpart at small $z$, but slightly larger at larger $z$ as they move away from region of high gravitational potential. For light rays emitted from the minimum density void, the differences between the GR and Newtonian calculations generally remain small—at the subpercent level—until the overdensity is approached. In the vicinity of the overdensity, the gravitational potential can be strong enough to cause a blue-shift, as evident in the top panel of Fig. 4.

Finally, we mention further details of the case with $\bar{\delta}=0.05$. This choice of initial conditions represents the extreme limiting case where the Newtonian treatment completely breaks down, and the Newtonian potential reaches $|\psi|\sim 1/2$ after a 15-fold increase of scale factor. As shown in Fig. 1, though the collapse at the overdensity (middle panel) occurs faster (in terms of proper observer time) in the Newtonian calculation than the full GR one, and the two calculations begin to noticeably differ well before halo formation, the evolution of the density in the void (top panel) still agrees well, with very little “backreaction” of the high-curvature region on the global expansion. In Fig. 5, we also show the luminosity distance-redshift relation for this case, which continues the trend found in Fig. 4, with increasing deviation between the Newtonian and GR calculations. Again, even for this extreme case, the differences between the light propagation in the void region are small. We are also not able to continue the GR calculation forward indefinitely, but it appears that a black hole is being formed at the overdensity. However, accurately tracking the attendant small scales requires adaptive mesh refinement, which we leave to future work.

IV Discussion and Conclusion

In this work, we have shown that a meaningful comparison can be carried out between standard $N$-body simulations of cosmological structure formation, which assume Newtonian-type gravity on the background of a homogeneously expanding universe, and full solutions of the Einstein-Vlasov equations, which make no assumptions regarding a background cosmology. For computational expediency, we have focused on a simple set of initial conditions, with inhomogeneities at a single length scale, but considered a range of amplitudes, including going all the way to the limit where the nonrelativistic assumptions underlying the Newtonian calculation break down. Tackling a more realistic power spectrum of density fluctuations will require more advanced techniques, such as adaptive mesh refinement, and will be quite computationally expensive given the stringent requirements place on time steps due to the fact that information propagates at the speed of light.

We find that for small initial density fluctuations, the Newtonian and GR calculations show excellent agreement (with differences typically subpercent and consistent with truncation error) well into the regime where the deviations from homogeneity become nonlinear. For large density fluctuations, the dominant relativistic correction seems to be that the collapse of overdensities occurs slower in the full GR calculation compared to the Newtonian one. These discrepancies can already be anticipated from the Newtonian calculation alone as the gravitational potential and infall velocities are approaching relativistic values. Even for such cases, the effect on the expansion outside the high density/velocity regions (e.g. in the voids) is found to be small, bounding backreaction effects.

Comparing the properties of light propagation in the Newtonian and GR calculations, we demonstrated that the resulting distance-redshift relations agree at the subpercent level as long as the Newtonian potential does not exceed the limit of a weak field approximation, i.e. $|\Psi_{N}|\leq 0.1$. As a limiting case, we have considered initial conditions all the way up to ones where the fluctuations in the density exceed the average value at the corresponding scales in the standard $\Lambda$CDM model by factor of $\sim 500$ (that is, at the present time they roughly correspond to $\sim 0.5$ at a Gpc scale). Since our simulations test the evolution on cosmological scales of perturbations with amplitudes exceeding those applicable to observational cosmology, we conclude that the obtained results provide a strong validation of the standard Newtonian approach employed in observational cosmology. In particular, our comparison implies that GR corrections to the Newtonian calculation of the cosmic variance in the local measurement of the Hubble constant are negligible. This strengthens the conclusion that a $\sim 9\%$ difference between the local and cosmic microwave background (CMB) based measurements of the Hubble constant, currently at $4.4\sigma$ statistical significance (Riess et al., 2019), cannot be ascribed to the cosmic variance which is estimated at $\sim 0.5$ percent (Wojtak et al., 2014; Wu and Huterer, 2017; Odderskov et al., 2014). This in line with the conclusion of a recent study in Ref. Macpherson et al. (2018) that looked at variations in the local expansion in a particular gauge using GR-fluid simulations (that hence cannot describe multistream regions) with a cosmologically motivated power spectrum.

The methods described here could be applied to study the formation of primordial black holes during a matter-dominated era (see e.g. Carr et al. (2017) and references therein), or scenarios where black holes make up some fraction of the dark matter. They could also be used to study ultralarge scale structure Braden et al. (2017), which could be related to understanding persistent CMB anomalies at large angular scales, which seem to indicate a violation of statistical isotropy and scale invariance of inflationary perturbations (Schwarz et al., 2016). Comparable scales will be also probed by the upcoming deep imaging cosmological surveys. In particular, the Large Synoptic Survey Telescope will reach an unprecedented effective volume of $\sim 4H_{0}^{-3}$ (LSST Science Collaboration et al., 2009).

Acknowledgements.

W.E.E. acknowledges support from an NSERC Discovery grant. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. R.W. was supported by a grant from VILLUM FONDEN (Project No. 16599). F.P. acknowledges support from NSF grant PHY-1912171, the Simons Foundation, and the Canadian Institute For Advanced Research (CIFAR). Computational resources were provided by XSEDE under grant TG-PHY100053 and the Perseus cluster at Princeton University.

Appendix: Numerical error results

In this appendix, we include some results on numerical convergence. For the GR simulations, we initially find the numerical error to be dominated by the grid spacing, which also sets the integration time step. However, at late times, as large under and overdensities develop, the number of particles used to sample the matter distribution becomes important. In Fig. 6, we show the convergence of the Einstein constraints with increasing numerical resolution for the $\bar{\delta}=10^{-2}$ case. The results have been scaled assuming second order convergence with grid spacing.

In Fig. 7, we show the halo crossing time for this same case as a function of resolution, for both the GR and Newtonian simulations. The discrepancies with resolution in the time of first crossing are small compared to the differences between the GR and Newtonian simulations (though they do become more pronounced for subsequent oscillations).

Finally, we compare the resolution dependence of the luminosity distance-redshift measures in Fig. 8. From this it can be seen that most of the $\lesssim 1\%$ differences between the GR and Newtonian simulations seen at early times or in the propagation outside the very high density regime are attributable just to truncation error. In contrast, the significant differences in propagation in the vicinity of the large overdensity exceed the truncation error, and in some cases are underestimated at lower resolutions.