Numerical Discreteness Errors in Multi-Species Cosmological N-body Simulations

TL;DR

This paper analyzes numerical discreteness errors in multi-species cosmological N-body simulations, identifying their sources and proposing mitigation strategies to improve the accuracy of large-scale power spectrum predictions.

Contribution

It provides a detailed analysis of discreteness errors in two-species cosmological simulations and introduces two effective mitigation methods to reduce these errors.

Findings

Discreteness errors cause biased small-scale power growth.

Large-scale power offsets persist even with equal particle masses.

Mitigation strategies significantly reduce power spectrum offsets.

Abstract

We present a detailed analysis of numerical discreteness errors in two-species, gravity-only, cosmological simulations using the density power spectrum as a diagnostic probe. In a simple setup where both species are initialized with the same total matter transfer function, biased growth of power forms on small scales when the solver force resolution is finer than the mean interparticle separation. The artificial bias is more severe when individual density and velocity transfer functions are applied. In particular, significant large-scale offsets in power are measured between simulations with conventional offset grid initial conditions when compared against converged high-resolution results where the force resolution scale is matched to the interparticle separation. These offsets persist even when the cosmology is chosen so that the two particle species have the same mass, indicating that the error is sourced from discreteness in the total matter field as opposed to unequal particle mass. We further investigate two mitigation strategies to address discreteness errors: the frozen potential method and softened interspecies short-range forces. The former evolves particles under the approximately "frozen" total matter potential in linear theory at early times, while the latter filters cross-species gravitational interactions on small scales in low density regions. By modeling closer to the continuum limit, both mitigation strategies demonstrate considerable reductions in large-scale power spectrum offsets.

Equations

$$r_h\approx 3\frac{L}{N_g}\equiv 3{\Delta }_G.$$

$${\Delta }_F=ϵ{\Delta }_P,$$

$${\Omega }_m$$

$$h$$

$$m_{\alpha }=\mathrm{bar}\rho {\Omega }_{\alpha }{\Delta }_P^3,$$

$${\delta }_{\alpha }(\text{k},z_i)=\mu (\text{k})\sqrt{P_{\alpha }(k,z_i)},$$

$$P_{\alpha }(k,z_i)=A{k}^{n_s}{T}_{\alpha }^2(k,z_i),$$

$$T_m^{\mathrm{bs}}(k,z_i)=\frac{D(z_i)}{D(z_f)}{T}_m(k,z_f),$$

$$T_{\alpha }^{\mathrm{bs}}(k,z_i)=\frac{T_{\alpha }(k,z_i)}{T_m(k,z_i)}{T}_m^{\mathrm{bs}}(k,z_i),$$

$$\text{x}(\text{q},z_i)=\text{q}+\text{d}(\text{q},z_i).$$

$${\text{d}}_{\alpha }(\text{k},z_i)=\frac{i\text{k}}{\Vert \text{k}{\Vert }^2}{\delta }_{\alpha }(\text{k},z_i).$$

$${\text{v}}_{\alpha }(\text{k},z_i)=\frac{i\text{k}}{\Vert \text{k}{\Vert }^2}{\frac{\partial {\delta }_{\alpha }(\text{k},z_i)}{\partial t}|}_{t=z_i}.$$

$${\text{v}}_{\alpha }(\text{k},z_i)=\frac{\dot{D}(z_i)}{D(z_i)}{\text{d}}_{\alpha }(\text{k},z_i).$$

$$T_{\theta ,\alpha }(k,z)\equiv -{\frac{d{T}_{\alpha }(k,z)}{dt}|}_{t=z}.$$

$$T_{\theta ,\alpha }^{\mathrm{bs}}(k,z_i)=\frac{H(z_i)}{a(z_i)}\frac{T_{\alpha }^{\mathrm{bs}}(k,z_i+\Delta z)-T_{\alpha }^{\mathrm{bs}}(k,z_i-\Delta z)}{2\Delta z},$$

$${\text{v}}_{\alpha }(\text{k},z_i)=-\frac{i\text{k}}{\Vert \text{k}{\Vert }^2}{\theta }_{\alpha }(\text{k},z_i),$$

$$\varphi (k,z_i)=-\frac{4\pi G\mathrm{bar}\rho }{k^2}{\delta }_m(k,z_i),$$

Full text

††pubyear: 2020††pagerange: Numerical Discreteness Errors in Multi-Species Cosmological N-body Simulations–12

1 Introduction

Numerical simulations of cosmological structure formation are an indispensable tool for providing theoretical predictions of the distribution of matter on non-linear scales. In combination with observational measurements, these predictions can be used to provide an understanding of structure formation and to solve the inverse problem of determining the cosmological makeup of the Universe. Next-generation observational probes – including DESI (DESI Collaboration et al., 2016), the Roman Space Telescope (Spergel et al., 2015), Rubin Observatory’s LSST (Ivezić et al., 2019), Euclid (Laureijs et al., 2011), and SPHEREx (Doré et al., 2014) – will continue to refine measurement constraints and thus require a matched increase in precision from theoretical modeling. From this point of view, there are two primary avenues that need to be advanced in simulations: 1) increased physical fidelity, including the accurate treatment of hydrodynamic forces and astrophysical processes impacting baryons, as well as inclusion of ingredients such as massive neutrinos, early dark energy, modified gravity, and primordial non-Gaussianity; 2) an improved understanding of the systematic errors inherent to the simulation approach, particularly with regard to the multi-species case that is both more numerically complex and less studied.

The earliest and most common form of cosmological simulation uses the N-body approach (e.g., Klypin & Shandarin, 1983; Efstathiou et al., 1985; Barnes & Hut, 1986; Hockney & Eastwood, 1988), which can be viewed as a dynamical Monte Carlo algorithm applied to the Vlasov-Poisson system of equations in an expanding universe, for which a set of point particles sampling the total matter field are evolved under their mutual gravitational potential. Each particle represents a parcel of the total matter fluid with a macroscopic mass that scales as the cube of the mean interparticle separation, $m_{p}\propto\Delta_{P}^{3}\equiv(L/N_{p})^{3}$, where $L$ is the simulation box width and $N_{p}$ is the number of particles along one dimension. It is common to solve the equations of motion using a gravitational force resolution, $\Delta_{F}$, that is finer than the mass resolution scale set by the mean interparticle separation, $\Delta_{P}$. In this case, it is natural to question whether the inherent discreteness of the particle representation leads to a numerical solution that systematically departs from the desired continuum limit. In the language familiar from plasma physics, we are interested in the circumstances under which the Klimontovich description (corresponding to the fine-grained distribution function) agrees with the mean field description underlying the Vlasov equation (Nicholson, 1983). We note that our concern here is not so much with particle collisionality effects (Okuda & Birdsall, 1970) which may be of concern in the deeply non-linear region, but with issues arising at much larger scales.

This question has received significant attention with various studies addressing the problem from a number of perspectives (e.g., Melott et al., 1997; Splinter et al., 1998; Knebe et al., 2000; Power et al., 2003; Diemand et al., 2004; Heitmann et al., 2005; Romeo et al., 2008; Joyce et al., 2009; Heitmann et al., 2010). The viewpoint normally taken is that discreteness errors are indeed introduced when $\Delta_{F}\lesssim\Delta_{P}$, but that these 1) remain confined to small scales due to the preferential transfer of power from large to small scales in the cold dark matter (CDM) paradigm, and 2) are effectively washed out at later times due to non-linear growth and mixing (Little et al., 1991; Hamana et al., 2002; Bagla & Prasad, 2009). Note that this is not necessarily the case in warm dark matter (WDM) or hot dark matter (HDM) scenarios where the lack of small-scale cosmological power is unable to mask discreteness errors, leading to obvious numerical artifacts that survive to late times (Götz & Sommer-Larsen, 2003; Wang & White, 2007; Melott, 2007; Power et al., 2016).

Discreteness errors become more readily apparent in particle-based simulations with multiple species as the symmetry amongst particles is broken and additional metrics can be analyzed. For instance, if the standard N-body method is modified so that the total matter field is split into two species of differing masses, discreteness errors manifest in the form of mass segregation within collapsed objects as energy is transferred from the heavy to the light species (Efstathiou & Eastwood, 1981; Binney & Knebe, 2002; Emberson et al., 2019). Such two-body interactions further extend to an artificial heating of gas (Steinmetz & White, 1997) in the situation where hydrodynamic forces are additionally modeled. In this case, one must also consider the more general interplay between the gravitational and hydrodynamic force resolutions to avoid numerical convergence issues, particularly when sub-resolution physics treatments are used (Bate & Burkert, 1997; Hopkins et al., 2018; Ludlow et al., 2019). In addition to modeling CDM and baryons, further discreteness effects occur when simulating neutrino particles. Recently, Sullivan et al. (2023) showed that the discrete signature of the initial condition neutrino particle grid can transfer onto the total matter field if the neutrino resolution is too coarse. Fortunately, in each of the above cases, the discreteness errors are mostly confined to small scales, making them easier to isolate and control.

Cosmological hydrodynamic simulations (e.g., Vogelsberger et al., 2014; Schaye et al., 2015; McCarthy et al., 2017; Emberson et al., 2019) tend to make the approximation that both the CDM and baryon tracer particles sample the total matter field (i.e., they are both initialized using the mass-weighted sum of the individual CDM and baryon transfer functions). This is technically incorrect since the disparate early universe evolution of the two species leads to markedly different transfer functions at the common initial redshifts ($z\approx 50-200$) of simulations. However, Angulo et al. (2013) showed that using the physically consistent procedure of initializing each species with its respective transfer function introduces a new discreteness error that, unlike the cases above, propagates to large scales. Specifically, if the particles are evolved with a force resolution $\Delta_{F}\lesssim\Delta_{P}$, then both the CDM and baryons acquire incorrect growth that extends all the way out to the largest scales at which the two transfer functions are initially different ($k\gtrsim$10^{-2}\text{,}\invMpch$$). The incorrect growth operates in opposite directions for the two species and has a greater amplitude in the lighter baryon particles. Angulo et al. (2013) attribute this problem to spurious particle coupling between the unequal mass species and suggest using an adaptive softening length in the CDM-baryon gravitational force interactions as a solution. Similar conclusions were also reached by O’Leary & McQuinn (2012).

Alternative solutions to this problem have been proposed by Bird et al. 2020 and Hahn et al. 2020, hereafter referred to as B20 and H20, respectively. The preferred method in B20 is to switch from the conventional offset grid pre-initial condition (pre-IC) setup to a mixed grid$+$glass (White, 1994) configuration. Another stated solution is to maintain an offset grid configuration, but to under-sample the baryon particles relative to the CDM so that each species has the same mass. In other words, the findings in B20 are consistent with the idea that the discreteness error arises from the fact that the particle masses are unequal. The method presented in H20 similarly involves modifications to the ICs. In this case, an offset grid configuration can still be used, but the individual transfer functions are encoded in mass perturbations instead of positional displacements. While the proposed methods in B20 and H20 successfully diminish the spurious growth of power, a complete theoretical explanation of both the propagation of discreteness errors and how the respective solutions formally mitigate these effects in a controlled way remains to be worked out.

In this paper, we continue the investigation into the origin of the large-scale discreteness error in mixed CDM plus baryon simulations. This is achieved using a controlled set of multi-species simulations that focus exclusively on gravitational forces so as to avoid the extra symmetry breaking that occurs when baryons are subject to hydrodynamic forces. We isolate for the impacts of having unequal masses versus unequal transfer functions by running simulations where either the mass or transfer function is held common for both species while the other is varied. We find that the discreteness error still exists in the case where individual transfer functions are used but the cosmology is chosen so that the particle masses are the same. In this case, the only difference compared to the standard cosmology is that the amplitude of the incorrect growth is shared equally between the CDM and baryon particles rather than mainly in the baryons. We argue that the discreteness error is not driven by having unequal particle masses, but more generally due to an unmixed discretization of the total matter field. In addition, this error is fundamentally distinct from artificial particle coupling that leads to small-scale features such as mass segregation.

In an effort to improve the discreteness error in multi-species cosmological simulations, we also evaluate various mitigation strategies. In addition to the proposed strategies in B20 and H20, as well as starting the simulation at later times, we propose two new strategies: 1) an approach based on evolving the particles under the linear total matter field at early times, and 2) using softened cross-species force interactions in low-density regions. We show that the biased power growth can be significantly reduced with either mitigation strategy.

The rest of the paper is organized as follows. We present an overview of our numerical methods and relevant background information in Section 2. In Section 3, we focus on simulations that use the standard approach of initializing the CDM and baryon components with the total matter transfer function and show that the discreteness error exists in this case, but is confined only to small scales. We then switch to the physically consistent initialization strategy in Section 4 where the two species are initialized with their respective transfer functions. Here we confirm previous findings that this leads to large-scale biases in the growth of both species and demonstrate that this occurs even if the particles have the same mass. This is further illustrated in Appendix A using a simple two-fluid plane wave toy model. In Section 5, we compare the proposed solutions from B20 and H20 and also present two new mitigation strategies to suppress the discretization error. Concluding remarks and further discussions are presented in Section 6.

2 Method

We provide below a summary of the numerical methods used in this work. We begin in Section 2.1 with a description of the cosmological simulations performed here, followed in Section 2.2 with further details on the multi-species ICs used to initialize each run. We finish in Section 2.3 with an outline of the frozen potential method that is later presented as a potential mitigation strategy for large-scale discreteness errors.

2.1 Simulations

The simulations in our paper are performed using the Hardware/Hybrid Accelerated Cosmology Code (HACC). HACC is a cosmological N-body solver designed to scale efficiently on all modern supercomputing platforms. We will briefly describe below the main code features of HACC relevant to this work and refer the reader to Habib et al. (2016) for more specific details. Despite the fact that we carry out mixed CDM plus baryon simulations, we restrict attention to purely gravitational forces so that we can easily identify particle asymmetries arising from discretization. In addition, we consider two different types of gravitational force solvers that are included in HACC to illustrate that discreteness errors arise regardless of the nature of the local force algorithm.

In the first method, which we refer to as HACC-PM, the gravitational force is computed using a spectral solver that employs a particle mesh (PM) algorithm to solve the Poisson equation in Fourier space. In this approach, particle masses are deposited to a uniform mesh containing $N_{g}^{3}$ cells via a cloud-in-cell (CIC) interpolation. In the second method, which we refer to as HACC-TPM, the PM code is used to evaluate a long-range gravitational force and is further augmented with a short-range force solver that traverses a tree data structure (Gafton & Rosswog, 2011) to compute direct pairwise gravitational interactions on small scales. In this case, a spectrally filtered “quiet” PM algorithm is matched to the short-range force at a handover scale, $r_{h}$, defined to equal roughly three grid separations:

[TABLE]

The short-range force is evaluated only for particle pairs confined to the handover scale and, for sufficiently distant leaves (determined via an opening angle criterion), a superparticle at the leaf center of mass is used instead of looping over all particles on the leaf. The short-range force is smoothed on small scales using a tunable Plummer softening to avoid artificially strong particle interactions. To expedite the particle search, we build a chaining mesh (CM; Hockney & Eastwood, 1988) containing cells of width $l_{\rm CM}=4\Delta_{G}$, chosen to encapsulate $r_{h}$, and construct an independent tree (“bushes”) within each cell.

We will consider simulations with a range of force resolutions in order to show that discreteness errors arise when the force resolution is made finer than the particle separation. In each run, we simulate $2\times N_{p}^{3}$ CDM plus baryon particles and define the mean interparticle separation of each species as $\Delta_{P}=L/N_{p}$. In the case of HACC-PM, the force resolution is set by the PM grid spacing, $\Delta_{G}$, which we adjust by varying $N_{g}$. For the HACC-TPM runs, we fix the PM mesh so that $N_{g}=N_{p}$ and instead vary the force resolution by setting the softening length, $r_{\rm soft}$, to be an adjustable fraction of the grid scale. In other words, the force resolutions of the two simulation methods can be expressed in terms of the mean interparticle separation as:

[TABLE]

where $\epsilon=N_{p}/N_{g}$ for HACC-PM and $\epsilon\approx r_{\rm soft}$ for HACC-TPM. We will show that discreteness errors arise whenever $\epsilon\lesssim 1$.

Throughout the paper, we assume a $\rm\Lambda CDM$ model with the following cosmological parameters:

[TABLE]

The corresponding mass of each simulation particle, $m_{\alpha}$, is set by the mean density of the Universe, $\mathrm{bar}{\rho}$, mean interparticle separation, $\Delta_{P}$, and energy density, $\Omega_{\alpha}$, of the given species:

[TABLE]

with $\alpha\in\{c,b\}$ denoting the CDM or baryon component, respectively. We consider two scenarios in our tests below: 1) an equal density cosmology with $\Omega_{c}=\Omega_{b}=0.15$ (yielding a particle mass ratio, $q\equiv m_{c}/m_{b}=1$, of the simulated CDM and baryon particles), and 2) $\Omega_{c}=0.25$ and $\Omega_{b}=0.05$, referred to as the “concordance cosmology” ($q=5$). The equal density case is used in Section 4.1 to investigate discretization errors in the absence of particle mass disparities. The concordance cosmology is used in all remaining runs and was designed to approximately match the energy density ratio observed with Planck (Planck Collaboration, 2020).

2.2 Initial Conditions

The general strategy for creating cosmological ICs is to start with a random realization of the target density field. This is achieved in Fourier space by taking the point-wise product of a Gaussian random white noise field, $\mu(\textbf{k})$, and the square root of the (dimensional) power spectrum:

[TABLE]

with $\alpha\in\{m,c,b\}$, where $m$ denotes the total matter field. Here $z_{i}$ is the initial redshift of the simulation which is taken to be $z_{i}=200$ for the majority of runs considered here. The power spectrum can be written in terms of the linear theory (LT) transfer function, $T_{\alpha}$, in the following manner:

[TABLE]

where $n_{s}$ is the spectral index and $A$ (not to be confused with the scalar amplitude $A_{s}$) is a normalization constant set by the value of $\sigma_{8}$.

There are generally two conventions in multi-species simulations for the choice of the target density field. The most common case is to use the total matter transfer function to initialize both the CDM and baryon particles with an implicit assumption that baryons follow the dark matter. We refer to this as the “single-fluid” case since both the CDM and baryon particles act as tracers of the same total matter fluid despite the fact that the two species are meant to individually represent the separate components of that field. The more formally consistent approach, especially at high redshifts, is to use the individual CDM and baryon transfer functions to initialize each respective species; we refer to this as the “two-fluid” case.

Transfer functions are normally computed using an Einstein-Boltzmann code (e.g., CAMB; Lewis et al., 2000). However, one needs to be careful of differences in the forward modeling assumed in cosmological simulations compared to that of Boltzmann codes. In particular, simulations implicitly work in the Newtonian framework of gravity, whereas Einstein-Boltzmann codes also include other physical effects such as the coupling between matter and radiation. Accordingly, the usual approach is to “back-scale” the transfer function evaluated at the final redshift of the simulation, $z_{f}$, to the starting redshift using the same forward modeling used in the simulation. In the single-fluid case, this process involves simple multiplicative factors:

[TABLE]

where $D(z)$ is the linear growth factor of the simulation forward model, normalized so that $D(0)=1$, and the superscript “bs” denotes the back-scaled result. This back-scaling strategy follows trivially from the fact that, in linear theory, the total matter density field evolves in a scale-independent manner, $\delta_{m}(\textbf{k},z)\propto D(z)$, in the Newtonian forward model.

This procedure needs to be modified in the two-fluid case since the density perturbations of CDM and baryons are both scale-dependent, even in the Newtonian forward model. For instance, baryons are coupled to photons via Compton scattering before recombination which suppresses the growth of the perturbations compared to CDM and also imprints baryon acoustic oscillations (BAO) on the baryon field. After recombination, the scale-dependent differences in the CDM and baryon density fields gradually transfer onto each other, though significant disparities still remain at the common starting redshifts used in simulations. The strategy used here is to take the logical extension of the standard back-scaling procedure in Eq. (7) one step further:

[TABLE]

where $\alpha=\{c,b\}$. In this way, we maintain the same total matter field as in the single-fluid case while ensuring that the relative proportions of the CDM and baryon density fields match the Einstein-Boltzmann predictions at $z_{i}$.

We use the back-scaled transfer functions, $T_{\alpha}^{\rm bs}$, to construct a random realization of the target density field following Eq. (5). This is evaluated on a uniform mesh containing $N_{p}$ cells along each axis. Next, we use the Zel’dovich approximation (ZA; Zel’dovich, 1970) to compute a displacement field, d, that maps particles from their Lagrangian location, q, to their starting position, x, in the simulation:

[TABLE]

The displacement field is derived from the density field in Fourier space (e.g., White, 2014):

[TABLE]

Unless otherwise specified, we use a pre-IC configuration where the CDM and baryon particles are placed on uniform grids that are offset from each other by $\Delta_{P}/2$. In this way, the two particle grids are maximally offset from each other; this configuration is chosen to minimize any artificial particle coupling at early times (Yoshida et al., 2003). In addition, we offset each particle grid by an amount $\Delta_{P}/4$ from the vertices of the IC mesh used to sample the displacement field. This is done to ensure that the two species are initialized in a symmetric manner with respect to the IC mesh (as well as the PM mesh used later in the force solver) so that the particle-grid relationships do not contribute to any of the particle asymmetries analyzed later. The displacement field is then spectrally interpolated at the particle positions (see Valkenburg & Villaescusa-Navarro, 2017). For glass configurations, we instead use an inverse-CIC interpolation in real space. In this case, we additionally deconvolve the displacement field with the CIC filter in order to compensate for interpolation smoothing.

The initial particle velocity is computed from the time derivative of Eq. (9) which, along with Eq. (10), yields

[TABLE]

In the single-fluid case, we can use the scaling $\delta_{m}(\textbf{k},z)\propto D(z)$ to show that

[TABLE]

Therefore, the velocity is in the same direction as the displacement field with a magnitude that scales with the time derivative of $D(z_{i})$. In other words, the initialization strategy of standard single-fluid simulations requires knowledge of only $T_{m}$, $D$, and $\dot{D}$; the latter two of which can be numerically evaluated for any given cosmology.

On the other hand, the two-fluid approach requires that we explicitly compute the scale-dependent time derivative of $\delta_{\alpha}$ in Eq. (11). In this case, it is convenient to define a velocity transfer function, $T_{\theta,\alpha}(k,z)$, where $\theta\equiv\nabla\cdot\textbf{v}$ denotes the velocity divergence field. The velocity transfer function follows from the linearized continuity equation (e.g., Peebles, 2021) as the negative of the time derivative of the density transfer function:

[TABLE]

In practice, the velocity transfer function can be computed using a finite difference of the synchronous gauge density transfer function sampled at $z\pm\Delta z$ (Valkenburg & Villaescusa-Navarro, 2017). For consistency, this must be performed using the back-scaled density transfer functions:

[TABLE]

where we choose $\Delta z=1.0$. The velocity transfer function is used to generate a random realization of the velocity field, $\theta_{\alpha}(\textbf{k},z_{i})$, following Eqs. (5) and (6), with the same white noise field and normalization constant that were used for the density field. The initial particle velocities become

[TABLE]

where the negative sign arises from the definition of the velocity transfer function.

We check for self-consistency in our two-fluid back-scaling approach by numerically integrating the initial CDM plus baryon density and velocity transfer functions defined in Eqs. (8) and (14) from $z_{i}=200$ to $z_{f}=0$ using the linearized Euler, continuity, and Poisson equations. The results are displayed in Fig. 1. In the top panel, we show the relative difference between the total matter transfer function inferred by HACC and that computed by CAMB at a variety of redshifts between $z_{i}$ and $z_{f}$. At early times, the simulated total matter density field will have systematically higher power on large sales. This difference gradually diminishes with decreasing redshift until the total matter transfer function matches CAMB on all scales at $z_{f}=0$. This is the intended result of the back-scaling procedure and matches the usual behavior exhibited in standard single-fluid simulations. In the bottom panel, we compare the ratio $T_{b}/T_{c}$ of the simulation to that predicted by CAMB. By construction, this ratio matches CAMB at $z_{i}$ and maintains excellent agreement at the $0.1\%$ level at all times for scales $k\lesssim$30\text{,}\invMpch$$. On smaller scales, we do observe a systematic departure at lower redshift which reflects the omission of Jeans damping in the simulation forward model. In any event, this discrepancy is rather small (at the $1\%$ level) and confined to scales that are heavily impacted by non-linearities and baryonic physics at late times.

2.3 Frozen Potential

We show later that discreteness errors manifest on large scales for two-fluid simulations run at high force resolution. The amplitude of these errors is greatest at high redshift when weak particle mixing amplifies the discrete nature of the simulated density field. An obvious way to reduce these errors is to start the simulation at lower redshift, but this is more difficult in the two-fluid case since no fully self-consistent higher-order strategy beyond the ZA currently exists (see H20 for a method that is second- or third-order in displacement and first-order in relative velocity). We will present an alternative approach based on the frozen potential (FP) approximation.

The FP method as an approximation to the full non-linear problem was introduced by Brainerd et al. (1993) and Bagla & Padmanabhan (1994) and is named after the fact that the linear gravitational potential, $\phi\propto\delta_{m}/a$, remains constant in an Einstein-de Sitter (EdS) universe where $\delta_{m}\propto a$. (For a further discussion and comparison with other approximate methods, see Vraštil & Habib 2020.) We use a modification to this scheme that accounts for the fact that the potential is not frozen in a general $\Lambda$CDM universe, but rather scales as $D/a$. The idea is that we first construct ICs in the usual manner described above with the additional step of storing in memory the initial gravitational potential,

[TABLE]

expressed here in Fourier space with comoving units. In this expression, $\delta_{m}(k,z_{i})$ is the total matter field of the given realization; i.e., evaluated with Eq. (5). Then, instead of evolving the particles under self-gravity, we scale $\phi(k,z_{i})$ by $D/a$ at each time step and use this to compute the gravitational acceleration, $\textbf{g}=-\nabla\phi$, evaluated on the PM mesh. In other words, we replace the discrete gravitational potential computed using the simulation particles with the continuous gravitational potential predicted by the linear total matter field. For consistency, short-range particle forces are omitted. This “frozen” potential method is used until some later handover redshift, $z_{h}$, chosen large enough that LT still mainly holds. Afterwards, the particles are evolved until $z_{f}$ using the usual force solver based on the particle potential field.

In Fig. 2 we compare LT power spectra to the output of a FP simulation evolved to $z_{h}=10.8$. The FP run is performed in a $512\text{\,}\Mpch$ box containing $2\times 256^{3}$ CDM plus baryon particles with an $N_{g}=1024$ PM grid. As can be seen in the top panel, the FP simulation accurately reproduces the LT power spectrum for both CDM and baryons on all scales down to roughly half the particle Nyquist frequency, $k\lesssim k_{\rm Nyquist}/2$, where $k_{\rm Nyquist}\equiv\pi/\Delta_{P}$. On smaller scales, the simulated power spectra systematically deviates from LT, owing to the finite particle resolution of the simulation. Prominent spikes are also seen on small scales, with the first at $2k_{\rm Nyquist}$, which simply trace the discrete Fourier modes of the IC particle grid. The total matter power spectrum ratio between the FP run and LT is shown in the bottom panel. All three curves overlap on large scales with fluctuations around the LT prediction. These fluctuations are imprinted from the random white noise field used to seed the ICs and denote the sample variance associated with the finite simulation volume. While the white noise profile of the total matter field remains relatively fixed during evolution (as demonstrated by the overlap of each curve in the lower panel), this is not quite the case for the individual CDM and baryon power spectra due to the scale-dependent mixing of the two species.

In Section 5.3, we treat the FP method as an alternative to using higher-order Lagrangian Perturbation Theory (LPT) to initialize two-fluid simulations at lower redshift. In this way, we start with the ZA (i.e., first order LPT) at $z_{i}$ and then evolve the particles under the FP method until $z_{h}$. We thus effectively extend the IC from $z_{i}$ to $z_{h}$, but in an improved manner over simply applying the ZA at $z_{h}$. The improvement comes from the fact that particle motion responds to the local potential at each timestep during the FP evolution, enabling particle trajectories to acquire curvature as opposed to remaining on straight lines with the ZA. As a result, the FP method better resolves shell crossing at later times and provides a more accurate treatment of the linear growth of small-scale structure (see Vraštil & Habib, 2020). The added utility of the FP method for the purposes discussed here is that it uses a continuous representation of the total matter field and therefore circumvents errors arising from the discrete particle field. The FP method is also computationally efficient since it only needs to evaluate the real-space gravitational force field, $\textbf{g}(z_{i})$, once and scales this result with $D/a$ at subsequent time steps. Of course, the FP approximation is based entirely on LT meaning that care must be taken on the choice of $z_{h}$. In addition, the force resolution is limited on small scales by the size of the PM mesh.

3 Single-fluid Simulations

In this section, we present the results of simulations where two species of unequal mass are initialized with the same total matter transfer function. This is the common approximation made in cosmological hydrodynamic simulations where separate CDM and baryon species are evolved from an initial realization of their combined field.

We begin in Section 3.1 with a demonstration that this setup leads to small-scale discreteness errors that manifest early during the simulation. The amplitude of this error increases with the mass ratio of the particles, but remains bounded in the asymptotic limit when the mass ratio approaches infinity (i.e., when one species becomes a passive tracer of the other species). We argue that this error arises from a coherent discretization of the total matter field into an unmixed representation of heavy and light particles. This is further corroborated in Section 3.2 where we use a Lagrangian refinement method (Hahn & Angulo, 2015) to show that the error is significantly reduced when the particles are interpolated below the force resolution scale.

3.1 Asymmetric Sampling

We present here results from a suite of 24 simulations that were initialized with two species on offset grids with displacements and velocities set by the total matter transfer function. We consider three simulation box sizes and resolutions: 1) $2\times 128^{3}$ particles in a $200\text{\,}\Mpch$ box; 2) $2\times 512^{3}$ particles in an $800\text{\,}\Mpch$ box; 3) $2\times 512^{3}$ particles in a $200\text{\,}\Mpch$ box. The first two configurations have the same mass resolution which is 64 times coarser than that of the third configuration. For each configuration, we perform three HACC-PM runs that differ in terms of the size of the PM grid with $N_{g}/N_{p}=1$, 2, and 4. In addition, we perform one HACC-TPM run with $N_{g}=N_{p}$ and a softening length $r_{\rm soft}=\Delta_{P}/8$. In other words, for each configuration, we run four simulations where the force resolution varies between $\Delta_{F}/\Delta_{P}=$ 1, 1/2, 1/4, and 1/8. Furthermore, for each force resolution, we perform one unequal-mass simulation with mass ratio $q=5$, as well as a reference simulation with uniform particle masses ($q=1$). In the latter case, we still use the same cosmological parameters (based on our concordance model) as the $q=5$ simulation, but simply set all particles to have the same mass. In what follows, we refer to the two particle species in the $q=5$ runs as “heavy” and “light” to emphasize that there is no physical distinction between the two groups given that they are initialized using the same transfer function and only gravitational forces are applied.

The results from these simulations are summarized in Fig. 3. The first three columns show the HACC-PM runs in ascending order of force resolution while the last column shows the HACC-TPM run. The rows denote the different simulation configurations, as indicated by the label in the final column. In each case, we show the power spectra of the heavy and light species from the $q=5$ run relative to the power spectra of the corresponding particle group from the equal-mass reference run. By symmetry, the power spectra of the two species in the equal-mass run are identical which means that comparing the $q=5$ and $q=1$ runs allows for a direct measure of how the evolution of each species changes when only the masses are varied, with all other parameters held fixed. Note that the absolute power spectra of each species depends on the force resolution of the given simulation, but our focus here is on examining only the relative change that occurs when the particle mass symmetry is broken in the $q=5$ runs.

When the force resolution is of the same order as the mean interparticle separation (i.e., $\Delta_{F}=\Delta_{P}$, as in the first column of Fig. 3), the unequal-mass power spectra remain in good agreement with the equal-mass results on all scales and redshifts. However, when the force resolution is made finer than the mean interparticle separation ($\Delta_{F}<\Delta_{P}$), the power spectra of the heavy and light species in the $q=5$ run systematically depart from their corresponding behavior in the $q=1$ run. In particular, in the two HACC-PM runs in the second and third columns, we see that the heavy (light) species exhibits a systematic enhancement (suppression) of power on small scales around the particle Nyquist frequency. This bias originates at very early times in the simulation, as demonstrated at $z=47$. The fact that this biased growth is clearly present when the particles are still rather separated from each other suggests that this effect is fundamentally different from particle scattering events such as mass segregation. The bias increases in amplitude as the simulation progresses, peaking around $z=10$ before subsequently falling in amplitude and being almost entirely erased by $z=0$.

Evidently, the qualitative trends are somewhat different in the HACC-TPM run where we instead see an initial suppression (enhancement) of the heavy (light) species for $k\lesssim k_{\rm Nyquist}$ which then reverses sign at $k\gtrsim k_{\rm Nyquist}$. The fact that differences arise between the two simulation strategies is perhaps not surprising given that their algorithms differ in their treatment of forces around the particle Nyquist frequency. Nevertheless, the main conclusion remains that the two species exhibit systematically different behavior that operates in opposite directions relative to the equal-mass reference run.

Comparing the first two rows of Fig. 3 allows us to check the dependence on simulation volume at fixed mass resolution. We find that the results are essentially identical in each column of the first two rows, indicating that the error depends only on the mass resolution of the simulation and is invariant to finite-volume effects. Conversely, comparing the top and bottom rows shows the dependence on mass resolution at fixed volume. At fixed force resolution, the main difference is that the discreteness error is almost entirely erased by $z=2$ in the bottom row compared to the top row which exhibits biased growth down to $z=0$. We attribute this to the finer mass and (absolute) force resolution in the bottom row that allows the simulation to resolve smaller scales and capture non-linear growth at earlier times. Non-linear growth aggressively mixes the particle species and seems to wash out the discreteness error that had accumulated earlier in the simulation. A similar argument can be applied to the trend seen when comparing the second and third columns in each row. The increased force resolution in the third column enhances particle mixing which reduces the overall amplitude of the discreteness error at any given redshift.

We can attempt to explain the origin of this discreteness error by considering the nature of the total matter density field seen by the force solver. In the unequal-mass simulations, we are effectively using the masses as a weighting function to up-sample (down-sample) the total matter field at the location of the heavy (light) species. As a result, the total matter field is not smooth below the particle scale and instead oscillates between peaks and troughs in a coherent fashion that follows from the regularity of the offset grids used in the IC. This coherence in the sampling strategy feeds back as a coherent pattern in the gravitational potential at the locations of the heavy and light particles. The result is an asymmetry in the gravitational force calculation that leads to biased growth operating in opposite directions for the two species. If the force resolution is made coarser than the particle scale, then the coherent peaks and troughs average out in the density field, precluding any biased growth. Alternatively, using pre-IC configurations that more adequately mix the two species (as explored in Section 5), one can diminish the discreteness error.

We can explore this idea further by considering the asymptotic limit of the tests above where the mass ratio is set to infinity ($q=\infty$). In this idealized scenario, the particle mass of one species is set to zero meaning that this species acts as a passive tracer of the gravitational potential sourced by the massive species. The result is shown in Fig. 4 where we use the first simulation configuration described above ($2\times 128^{3}$ particles in a $200\text{\,}\Mpch$ box) performed using HACC-PM with a PM mesh so that $\Delta_{F}=\Delta_{P}/2$. As before, we normalize the results to the equal-mass reference simulation evolved under identical settings. For comparison, we also display the $q=5$ run shown previously in the second column of the top row in Fig. 3.

We find qualitatively similar results to the $q=5$ simulations. Namely, we see that the passive species displays heavily suppressed growth on small scales relative to the massive species. We can attribute this to the systematic offset of $\Delta_{P}/2$ between the locations at which the massive species sources the density field and the locations at which the passive species samples the gradient of the gravitational potential. By symmetry, we would instead find identical growth between the two species if the passive particles were initialized on top of the massive particles. Comparing to the $q=5$ case shows how this scenario changes when we transfer some of the mass from the massive to the passive species. The result is that the light particles still exhibit suppressed growth relative to the heavy species, but with a reduced amplitude compared to the $q=\infty$ case. Finally, symmetry dictates that the scenario changes in the $q=1$ case with both species growing in an identical manner.

To summarize, these tests demonstrate that using unequal particle masses to coherently up- and down-sample the total matter field leads to a corresponding growth bias operating in opposing directions for the two species. Relative to the symmetric equal-mass case, we find that the heavy (light) particles acquire enhanced (suppressed) growth on scales around the Nyquist frequency when evolved with a PM solver. The situation changes slightly when pairwise forces are used on small scales, but the main conclusion of asymmetric biased growth remains. This error stems from the fact that the force solver sees a total matter density field that is coherently discretized into unmixed components. Eventually, the error stops growing when the particle species become sufficiently mixed and the coherence in the density field is smoothed out. Subsequent non-linear growth tends to overtake the initial accumulation of this error.

3.2 Lagrangian Refinement

To further illustrate that the biased growth arises from an unmixed sampling of the total matter density field, we utilize the Lagrangian refinement method of Hahn & Angulo (2015) to smooth the particle discretization for a fixed force resolution. Each pre-IC grid cell is effectively replaced by $4^{3}=64$ refined particles uniformly distributed within the Lagrangian cube. We perform this upsampling on-the-fly during each force calculation step, without having to run the entire simulation at a higher particle resolution. This is achieved through a continuous mapping between Lagrangian space and the phase-space distribution, $\textbf{q}\mapsto(\textbf{x}_{q}(t),\textbf{v}_{q}(t)$), where we utilize a tri-linear interpolation map during the CIC deposition step in our integrator. While we do not expect the refining procedure to resolve the numerical solution at higher fidelity, we do expect it to mitigate discretization errors as the density sampling is now being smoothed.

In Fig. 5 we compare power spectra from one of the $q=5$ simulations performed both with and without the Lagrangian refinement method. Each simulation evolved $2\times 128^{3}$ particles in a $200\text{\,}\Mpch$ box using HACC-PM with $\Delta_{F}=\Delta_{P}/2$. In the top panel, we present the power spectrum of the heavy and light species from the two $q=5$ runs relative to the $q=1$ reference simulation; the bottom panel compares the total matter power spectra. We see that the Lagrangian refinement significantly reduces the level of biased growth in both the heavy and light species. In contrast, the total matter power spectrum is heavily suppressed relative to the reference run.

The Lagrangian refinement run differs in only the fact that the particle density field is interpolated to a scale 4 times finer than the mean interparticle separation. The result is a significant reduction in the biased growth on small scales that is achieved at fixed force resolution without increasing the actual number of simulation particles. This result reinforces the notion that the biased growth is attributed to the force solver resolving an unmixed particle distribution. We note that the Lagrangian refinement method is accurate only as long as the phase-space submanifold is well-sampled by the $N$-body particles (e.g., at early times in WDM simulations, see Hahn & Angulo, 2015). Once the particle sampling is no longer sufficient to trace the folding of the dark matter sheet in phase-space, the now unreliable mapping from Lagrangian to Eulerian space leads to a suppression of growth near the Nyquist frequency, as shown in the bottom panel of Fig. 5. Given this drawback, we do not recommended using Lagrangian refinement to address discreteness errors, but rather present it here as an illustrative example of how the biased growth is fundamentally sourced from the discrete nature of the density field in multi-species simulations.

4 Two-fluid Simulations

We now investigate the consequences of initializing simulations with individual density and velocity transfer functions for CDM and baryon particles. This is the more physically consistent approach for cosmological simulations despite being rarely used in practice. In order to separately distinguish numerical artifacts seeded by particle mass discrepancies (as investigated in Section 3), we begin by studying uniform mass ($q=1$) simulations in Section 4.1. In this case, we find that both species exhibit biased growth on large scales reinforcing the notion that the error is sourced from the unmixed discretization of the total matter field into its CDM and baryon components as opposed to arising purely from having unequal particle masses. We then further study unequal-mass ($q=5$) simulations in Section 4.2, which are found to behave in a similar fashion as the equal-mass case, but with the discreteness error shifted mostly into the lighter baryon component.

4.1 Balanced Cosmology

We begin with an analysis of a set of simulations that employ an equal density cosmology with $\Omega_{c}=\Omega_{b}$ so that the resulting CDM and baryon particle masses are equal. A total of four simulations were run, with each evolving $2\times 256^{3}$ particles in a $512\text{\,}\Mpch$ box. As before, we performed three HACC-PM runs that vary the PM mesh size so that $\Delta_{F}/\Delta_{P}=1$, 1/2, and 1/4, as well as one HACC-TPM run where the PM mesh is fixed with $N_{g}=N_{p}$ and the softening length is set so that $\Delta_{F}/\Delta_{P}=1/8$. The results are shown in Fig. 6 with each column denoting an individual run while the rows separately show the baryon, CDM, and total matter power spectra, respectively.

The power spectra of each simulation is shown relative to a reference run evolved under the FP approximation. The reference run was performed using HACC-PM with $\Delta_{F}=\Delta_{P}$ and, as described in Section 2, evolved particles under the LT total matter field at each time step. We used the FP run as a baseline as opposed to simply comparing against LT since we found that the white noise profiles of the individual CDM and baryon power spectra become mixed when evolved together, obfuscating the signatures that we are looking for on large scales. While the FP method achieves the correct linear growth of each species on large scales (see Fig. 2), it does omit non-linear terms meaning that it is a reliable baseline only at high redshift. In addition, the FP method differs from the actual simulations in that it does not use any interpolation function to compute the density field. As such, we cannot make an apples-to-apples comparison between the simulations and the FP baseline on small scales. Neither of these issues are important for our application, however, since we are mainly focused on large scales and early times.

The results in Fig. 6 are analogous to those found earlier for the single-fluid simulations in Fig. 3. In particular, we find that the CDM, baryon, and total matter power spectra match the growth of the FP method on large scales for the runs in the leftmost two columns where $\Delta_{F}/\Delta_{P}\geq 1/2$. Conversely, we observe a striking contrast in the growth of the CDM and baryon components in the rightmost two columns for which $\Delta_{F}/\Delta_{P}\leq 1/4$. Here we see that the CDM (baryons) grow systematically fast (slow) relative to the baseline run. Unlike the single-fluid simulations, the biased growth in the two species is clearly visible out to the largest scales resolved in the simulation. Moreover, we see prominent wiggles in both the CDM and baryon curves on large scales. This reflects a modulation in the BAO signal that results from a delayed transfer of BAO power from baryons to CDM due to the biased growth of each species. Given that the errors extend to large scales, we do not observe any significant differences between the HACC-PM and HACC-TPM runs which differ only in their treatment of small-scale forces. Interestingly, we find that the biases of the individual components in each simulation cancel out in the total matter field which is seen to match the large-scale growth of the reference run.

We argue that the errors seen in the single- and two-fluid simulations are manifestations of the same problem; namely, the unmixed discretization of the total matter field. In the single-fluid case, the confinement of the error around the Nyquist mode reflects the fact that the discretization occurs only at the particle scale. In the two-fluid case, the error extends to the much larger scales where the individual transfer functions differ from that of the total matter. Indeed, in simulations with much larger boxes, we have found that the bias does not materialize on scales where the CDM and baryon transfer functions converge ($k\lesssim$10^{-2}\text{,}\invMpch$$). Importantly, the bias occurs even in simulations with uniform particle masses. This confirms that the large-scale error seen in two-fluid simulations is fundamentally sourced by the discretization of the total matter field and not from particle masses per se. We reinforce this interpretation using a plane wave case-study in Appendix A, exploring discreteness effects in an idealized setting.

4.2 Concordance Cosmology

The previous subsection served as an illustrative example of how large-scale discreteness errors still occur in simulations with uniform particle mass. We switch attention here to the more realistic concordance cosmology which leads to unequal-mass ($q=5$) CDM and baryon particles. This is similar to the mass ratio used in previous works (Angulo et al. 2013; B20; H20) that have also studied discreteness errors in two-fluid simulations.

For this test, we ran a single HACC-PM simulation that evolved $2\times 256^{3}$ particles in a $512\text{\,}\Mpch$ box with $\Delta_{F}=\Delta_{P}/4$. We also performed a baseline run where the particles are evolved under the FP approximation with $\Delta_{F}=\Delta_{P}$. We compare the power spectra from these two runs in Fig. 7. This plot can also be compared to the third column of Fig. 6 where the same set of parameters (apart from the cosmology) are used. We observe here the same qualitative features as seen previously for the equal density cosmology. More specifically, the large-scale offsets and BAO modulation are present in both the baryon and CDM power spectra. The main difference, however, is that the bias in the CDM power spectra is much more subdued in the concordance cosmology. We attribute this to the fact that the CDM now dominates the total matter field meaning that the error is preferentially shifted into the subdominant baryon component. We again find that the large-scale growth of the total matter field aligns well with the reference run.

These results are consistent with the notion that errors arise from the unmixed discretization of the total matter field. In this case, the total matter field is split into particles with asymmetries in both mass and power spectra. The results of Section 3 suggest that the particle masses impart errors on small scales while the findings in Section 4.1 confirm that the disparate density fields are responsible for driving errors to large scales. In the next section, we examine various methodologies that can be used to address these issues.

5 Evaluation of Mitigation Strategies

We conclude our analysis with an evaluation of several approaches to address the large-scale discreteness errors seen in two-fluid simulations. For the remainder of this section, we utilize target simulations evolving $2\times 256^{3}$ particles in $64\text{\,}\Mpch$ boxes. This setup was chosen so that the mitigation strategies are evaluated using a mass resolution that is more representative of production-scale simulations used to generate synthetic skies. With our concordance cosmology, the CDM and baryon particle masses are $m_{c}=$1.1\times 10^{9}\text{,}\massh and $m_{b}=$2.2\times 10^{8}\text{\,}\massh, respectively. The simulations are performed using either HACC-PM with an $N_{g}=1024$ PM mesh ($\Delta_{F}=\Delta_{P}/4$) or HACC-TPM with $r_{\rm soft}=\Delta_{P}/20$. We have also performed a reference HACC-PM run that was chosen to have the same force resolution as the HACC-PM target runs but containing $2\times 1024^{3}$ particles ($\Delta_{F}=\Delta_{P}$) so that large-scale discreteness errors are absent. The reference run was initialized using the same white noise field as the target runs.

In Section 5.1, we investigate two initial condition approaches recently proposed by B20 and H20, and evaluate how well they reduce the offset power bias. We further analyze simulations evolved at lower starting redshifts in Section 5.2, where the sampling distribution is more mixed. Lastly, in Sections 5.3 and 5.4 we discuss new methods that can mitigate the discreteness errors while using standard grid ICs.

5.1 Initial Condition Configurations

In this section, we investigate the effects on the power spectrum when we start simulations with different ICs. Three separate initialization layouts at $z_{i}=200$ are considered here as follows: offset grid, mixed grid plus glass (B20), and mass perturbed (H20) ICs. Briefly, the configuration proposed by B20 is designed to efficiently mix the two species by initializing the CDM particles on a grid while the baryons use a glass distribution (White, 1994). The method described in H20 utilizes offset grids by displacing particles with the total matter transfer function, and therefore, avoids issues with unmixed discretization of separate fluids (with the additional benefit of optionally initializing to higher order LPT). The individual species transfer functions are separately encoded using mass and velocity perturbations. We use the publicly available MP-GenIC (Feng et al., 2018) code to generate the B20 ICs and monofonIC music-2 111Available from https://bitbucket.org/ohahn/monofonic for the H20 runs.

The results of three HACC-PM simulations initialized with the methods described above are presented in Fig. 8. From top to bottom, the panels show the power spectrum ratios between the simulations and the reference run for baryons, CDM, and the total matter. For the offset grid ICs, we find the previously discussed deviations (see Section 4) in the baryon and CDM power spectrum at all $k$, peaking at $z\sim 10$ with approximately $10\%$ suppression for baryons and $2\%$ enhancement for CDM. At later times, the large-scale offsets decrease as non-linear growth starts to govern the evolution at high $k$. Note that the discreteness errors only affect the growth of the individual species; the total matter power agrees with the reference simulation at low $k$.

In comparison, the grid$+$glass runs show significantly improved power spectra results, where the large-scale differences are within $2\%$ for baryons and less than $1\%$ for CDM at all redshifts. This reflects the pseudo-random nature of the glass which prevents coherence in the discretization of the total matter field into its two components. The drawback of the glass, however, is the introduction of noise on small scales, manifesting as an excess of power at medium and high $k$ in the individual and total matter distributions. We note that B20 has also shown that undersampling the baryonic fluid to maintain equal mass between the species results in similar improvements to a mixed grid$+$glass IC. We have confirmed this to be true using an offset grid simulation containing $256^{3}$ CDM and $150^{3}$ baryon particles. In accordance with our previous results, this improvement likely arises from symmetry breaking due to the irregular grid separation between species and not from having the same particle masses. Given that equal particle sampling of both fluids is the most common method of discretization, we focused on that configuration throughout this study.

Turning to the results of the mass-perturbed simulations, we find that the power spectra displays roughly $3\%$ suppression for baryons at $z=151.3$ and less than $1\%$ enhancement for CDM at all four redshifts. The baryon suppression quickly drops to the sub-percent level for $z\leq 10$. The success of the mass-perturbed method can be attributed to the fact that the individual transfer functions are absorbed into the masses as opposed to the initial displacements. In this way, the force solver sees a particle distribution that is consistent with the total matter field on large scales, but with small mass perturbations at the particle scale. Over time, the separate initial velocity fields will lead the particles to “dephase” meaning that the original issue causing the large scale offsets will resurface (i.e., see the small bias in baryon power at $z=10.8$ that is later reduced by non-linear growth). For this reason, the mass-perturbed method is more aptly suited for lower redshift initializations, as we show in the next section.

5.2 Low Redshift Initialization

A direct approach to address errors incurred by evolving unmixed distributions is achieved by initializing the simulation at lower redshifts, perturbing the particles further from their initial grid positions. A major concern about late simulation starts is the requirement of higher order LPT to accurately capture the curved trajectories of the initial displacement. In the single-fluid case, studies have illustrated at least second order perturbation theory (2LPT) is desired for initializing accurate N-body simulations at $z_{i}=50$ (e.g. Schneider et al., 2016; Michaux et al., 2020). However, there is not yet a consistent method to compute higher order LPT for two-fluid simulations. The procedure outlined in H20 does utilize higher order LPT for the combined baryon and CDM distribution, albeit using an initial relative velocity field, $v_{\rm bc}$, calculated only to first order. We will use both mass-perturbed (initialized with first and second order LPT) and standard grid configurations to study low redshift initial conditions below.

In Fig. 9, we show the power spectrum results from HACC-PM simulations initialized with three different methods at $z_{i}=46.9$: a standard grid setup with first order LPT and mass-perturbed runs with both first and second order LPT. Starting with the standard grid runs, we notice that although the large-scale offsets in the baryon and CDM power spectrum are still present, they are significantly reduced compared to the run initialized at $z_{i}=200$ (Fig. 8). More specifically, the power suppression for baryons on large scales lowers from approximately 10% to 3% at $z=10.8$. Next, the simulation initialized using the mass-perturbed method of H20 with 1LPT further improves the large-scale behavior of the individual species power, matching the high resolution reference simulation at the lowest $k$ modes. While the individual power spectra for the 1LPT standard grid and mass-perturbed runs differ on large scales, their total matter power remain in broadly good agreement on the full set of scales and redshifts. On mid to small scales, however, there is a clear lack of total matter power compared to the reference simulation. This discrepancy is significantly improved using the mass-perturbed 2LPT initial condition which displays sub-percent agreement with the reference run at $z=0.5$. Hence, the mass-perturbed 2LPT run not only removes the large-scale offsets, but also more closely matches the reference simulation at intermediate scales.

Thus, we have confirmed that lower redshift starts can reduce the biased growth in power on large scales. Moreover, the mass perturbation approach further addresses the necessity for higher order LPT initial conditions, in an improved manner over simple grid runs. Of course, if one wishes to study high-redshift structure formation (i.e., first galaxy formation) for which early start simulations are needed, consideration of the necessary particle sampling to overcome these effects is required. We exclusively utilize power spectrum diagnostics in this paper but there are in principle other metrics (e.g., bispectrum, halo mass function) that could be impacted by initializing the simulation at lower redshifts. Further studies regarding this topic are currently in progress.

5.3 Frozen Potential Approach

We present here an alternative method for starting offset grid simulations at low redshift while still maintaining a first-order initialization. The idea is to construct 1LPT initial conditions in the usual manner at $z_{i}=200$, followed with subsequent evolution to some later handover redshift, $z_{h}$, under the FP approximation (see Section 2.3). The remaining evolution is completed with the full force solver. In this way, the FP method can be thought of as an extension to the initial conditions. Of course, care must be taken on the choice of $z_{h}$ since a value too high will still suffer from discreteness errors while a value too low will delay non-linear growth. In general, the compromise between these limits will depend on the resolution of the simulation.

To explore this further, we perform three runs with our target simulation where the FP method is used until $z_{h}=46.9$, $21.3$, and $10.8$. In each case, HACC-PM is used to evolve the particles from $z_{h}$ to $z=0$. The results of these tests are organized in Fig. 10. The simulation with the highest handover redshift, $z_{h}=46.9$, shows approximately $3\%$ large-scale suppression in baryon power spectrum at $z=10.8$ and no significant large-scale offsets in CDM and total matter power spectrum at any redshift. This is a significant improvement compared to the $10\%$ suppression for baryons and $2\%$ enhancement for CDM shown in Fig. 8 for the default target run with $z_{i}=200$. On intermediate scales, however, there is a clear suppression in both baryon and CDM power for $z\geq 2$. This power suppression is replaced at $z=0$ with random noise at the few-percent level. At this point, the simulation volume is almost entirely non-linear which appears to mostly erase the earlier systematic power suppression. Lowering the handover redshift to $z_{h}=21.3$, the large-scale deviation in the baryon power spectrum further reduces to less than $2\%$. Unsurprisingly, the power suppression on intermediate scales has become more pronounced, though this is once again mostly erased by $z=0$. Finally, the simulation with $z_{h}=10.8$ shows significant power suppression at all times indicating that such a low handover redshift is not appropriate for capturing non-linear growth.

In summary, we find that using the FP approximation to evolve particles from first-order ICs at high redshift to later times can significantly reduce large-scale biases resulting from insufficient particle mixing. The challenge is that the FP method is limited by both its assumption of linear growth as well as a coarser force resolution due to the omission of short-range forces. In general, a careful analysis of the appropriate choice of $z_{h}$ would need to be performed for a given mass resolution, in addition to consideration of other metrics outside of the power spectrum.

5.4 Softening Interspecies Interactions

Our final mitigation strategy for reducing the large-scale power bias is connected to the observation that these errors occur when the force resolution is finer than the interparticle separation. In particular, errors are most pronounced in undersampled regions where the unmixed initial grid configuration of each species maximally impacts the evolution. In response, we explore a modification to our short-range gravity solver which adaptively turns off or softens interspecies interactions in low-density regions.

Within HACC, the short-range solver operates on neighboring chaining mesh (CM) cells of size $l_{\rm CM}=4\Delta_{G}$ (see Section 2.1). Utilizing the same spatial decomposition, we implement a particle number per cell threshold, $N_{t}$, below which the short-range cross-species interaction forces are turned off or softened for any particle within the cell. These include pairwise force contributions from cell neighbors. This modification only affects low density regions within the simulation, with the intention to not suppress non-linear growth in clustered regions. Note that for locations of the simulation at average background density we expect $N_{\rm CM}=2\times 4^{3}=128$ particles per cell. We found for our simulation resolution that a value of $N_{t}=160$, corresponding to an overdensity $\delta_{\mathrm{CM}}=0.25$, to be a good compromise between alleviating the discreteness errors while preserving the growth of non-linear structures and we adopt this threshold in the following.

We implement three different variations of the short-range force operator: SR-SOFTEN, SR-OFF, and SR-BOOST. SR-SOFTEN is motivated by previous mitigation attempts using adaptive gravitational softening lengths (e.g., Angulo et al., 2013; Villaescusa-Navarro et al., 2018) and the observation that the discreteness errors disappear when using HACC-PM with $\Delta_{F}=\Delta_{P}$ (see Section 4). In cells below the particle threshold $N_{\rm CM}<N_{t}$, we increase the softening length from $r_{\rm soft}=\Delta_{P}/20$ to $r_{\rm soft}=\Delta_{P}$ for cross-species interactions within the cell and its neighbors, effectively reverting the gravity solver to HACC-PM with $\Delta_{F}=\Delta_{P}$ for those interactions. SR-OFF extends that modification by completely turning off cross-species short-range forces, thereby isolating the effects of interspecies interactions on structure formation. Finally, in an attempt to compensate for the suppression of growth on small scales, we implement SR-BOOST, which boosts the particle masses by an amount $\Omega_{m}/\Omega_{\alpha}$ in same-species interactions when cross-species forces are turned off to match the total matter density. Note that the PM force calculation is unmodified in each of these cases meaning that cross-species interactions are still present on large scales.

The power spectrum ratios between the simulations with modified interspecies forces and the high-resolution HACC-PM reference simulation are shown in Fig. 11. We also include the HACC-TPM simulation with constant softening length ($r_{\rm soft}=\Delta_{P}/20$) and the HACC-PM simulation with $\Delta_{F}=\Delta_{P}$ for comparison. As previously discussed, the HACC-TPM simulation displays significant offsets of the large-scale CDM and baryon power at higher redshift (i.e., $z=10.8$) and extra power on small scales at lower redshift (i.e., $z=2$ and $z=0$) compared to the reference simulation, whereas the HACC-PM simulation does not show any offsets on large scales but has suppressed power on small scales because of the coarse force resolution. The SR-SOFTEN variation significantly improves the large-scale power spectrum for both species: while CDM no longer shows any notable large-scale offset, the baryon power is still about 1% suppressed at $z=10.8$, compared to the 12% suppression in the HACC-TPM run. When we completely remove the short-range cross-species interaction in the SR-OFF option, the offsets on large scales disappear for both CDM and baryons, confirming our initial statement that the large-scale offsets likely arise from the cross-species communication. In both the SR-SOFTEN and SR-OFF simulations, small-scale structures initially grow at a lower rate than in the HACC-TPM simulation, resulting in less power on large $k$ at higher redshift (e.g., $z=10.8$). This is due to the over-smoothing (SR-SOFTEN) and the removing (SR-OFF) of a subset of the short-range forces, essentially reducing the amplitude of the local density perturbation each particle experiences.

The SR-BOOST variant is added in an effort to compensate the lost power from the missing cross-species interactions by boosting the “remaining” short-range forces to the total matter density expected if baryons would perfectly trace CDM. Fig. 11 shows that the SR-BOOST method retains the low-$k$ results of the SR-OFF variant and does not suffer from discreteness errors. On medium and small scales, the power spectrum is significantly improved. In particular, we can directly compare the total matter power spectrum to the HACC-TPM simulation, for which $P_{m}$ is not affected by the large-scale discreteness errors. We find only small differences at $z=10.8$ and good agreement at $z=2$ and $z=0$.

We conclude by emphasizing that the SR-BOOST methodology assumes an idealized approximation that both fluids identically trace the total matter density. This picture is not only inconsistent in multi-fluid simulations, but also is further violated when hydrodynamic forces are additionally included. The SR-BOOST results are solely intended to highlight the cause of the discreteness errors, where directly removing cross-species interactions can drastically improve the power biases in the simplified gravity-only case we explore here. On the other hand, adaptive force softening methods – akin to our SR-SOFTEN implementation – are more generally applicable to multi-species simulations (and hence commonly used in cosmology codes), and are useful mitigation techniques to reduce the erroneous growth, as our results confirm.

6 Conclusion

The statistical precision and accuracy of survey observations are increasing steadily, challenging the required predictive capability of cosmological simulations. As the demands on error control become tighter, the long-standing issue of generating accurate multi-species initial conditions needs to be examined closely, especially given the low initial particle sampling density characteristic of cosmological simulations. In this paper, we investigate numerical discreteness errors in gravity-only simulations, where interspecies fluid ICs and evolution can be studied in a simplified setting. Our analysis used measurements of the density power spectrum, a commonly employed diagnostic particularly suited to studies of structure formation irregularities.

We have found that coherently unmixed discretizations of the total density field – encountered when initializing fluids with simulation tracer particles on offset grids – lead to biased growth in the individually sampled components if the force resolution employed is finer than the particle spatial sampling scale, i.e., $\Delta_{F}<\Delta_{P}$. In single-fluid simulations where the total field is split into species with non-uniform mass, we found that the error is confined to small scales around the (particle sampling) Nyquist mode. This behavior can be simply interpreted by examining the asymptotic limit where one species is massless; the suppressed growth of the tracer fluid is understandable given the systematic positional offset from the massive species particles that source the gravitational force. In two-fluid simulations, biased growth occurs out to the largest scales at which the individual transfer functions of the species differ, regardless of whether or not the particle masses of the two species are non-uniform. In both the single and two-fluid cases, the error arises immediately in the simulation and continues growing until the particles become sufficiently mixed due to non-linear gravitational evolution. Subsequent non-linear growth may reduce this error, but significant large-scale biases may still exist at $z=0$ if proper mitigation strategies are not implemented.

We evaluated several procedures for reducing discreteness errors in two-fluid simulations. First, we examined two previously proposed strategies that involve modifications to the initial condition setup. We found that the large-scale biases in the power spectra for both CDM and baryons are significantly reduced by starting with mixed grid plus glass (B20) or mass-perturbed (H20) ICs. The pseudo-random glass distribution successfully destroys coherence in the total matter discretization while the mass-perturbed method shifts the discretization to small scales at early times. Additionally, it is generally advantageous to initialize at later redshifts when the two species are better mixed. We showed that the biased growth on large scales is markedly improved with $z_{i}\sim 50$ starts for both the standard grid and mass-perturbed ICs. In this case, it is important to use higher order LPT, though multi-species methods incorporating relative velocities beyond first order do not currently exist.

We considered a frozen potential approach as an alternative to starting at lower redshifts while still maintaining offset grid initial conditions. This methodology enforces linear theory evolution up to a specified time, and was also successful in reducing the measured discretization effects. However, avoiding significant suppression of non-linear growth as well as specifying a resolution-dependent redshift parameter complicate its usage. A final mitigation strategy we examined modified the calculation of short-range interspecies gravitational forces in low-density regions. For example, we found that turning off interspecies interactions within chaining mesh cells below a threshold density effectively removed the large-scale offsets for both CDM and baryons. In general, adaptive softening of gravitational forces reduces the discreteness errors and is a viable approach for addressing this effect.

The focus of our investigation was entirely devoted to power spectrum analyses, and therefore does not assess the impact of multi-species numerical discreteness errors on other important metrics such as the halo mass function or bispectrum. More work is needed in these areas in order to obtain a deeper understanding of the full impact of this simulation artifact. In addition, the simulations presented here were restricted to two species evolved purely under gravity (and, in most cases, treated symmetrically by the force solver). Discreteness errors have the potential to become more pervasive if coupled to a hydrodynamic solver or if a third species such as neutrinos are included. Furthermore, the impact of these effects are not clear when other types of multi-species solvers are considered (e.g., Eulerian or hybrid methodologies). In order for cosmological simulations to meet the stringent prediction demands of near-future observational surveys, it is of utmost importance to fully understand and better control the systematics associated with the discrete nature of particle-based simulations.

Acknowledgments

We thank Oliver Hahn, Katrin Heitmann, and Patricia Larsen for helpful discussions and comments on the draft. SH additionally thanks Adrian Melott, Robert Ryne, Sergei Shandarin, and Martin White for past discussions on errors due to discreteness effects. Argonne National Laboratory’s work was supported under the U.S. Department of Energy contract DE-AC02-06CH11357. This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award cusp in 2021/22.

Data Availability

The analysis data from the simulations presented in this article will be made available upon request.

Appendix A Two-fluid Plane Wave Experiment

We provide here an illustrative example of how an unmixed discretization of the total matter field leads to biased growth in two-fluid simulations. This is accomplished using an idealized setup where two particle types (“CDM” and “baryon”) of uniform mass are initialized using a single plane wave density perturbation with a species-dependent amplitude. We find that the species with the larger (smaller) perturbation amplitude exhibits enhanced (suppressed) growth due to paired particle coupling producing an inward (outward) force around peaks in the total matter density field. If we add perturbations to other scales, but with an equal amplitude for both species, then the biased growth remains confined only to the scale with the species-dependent amplitude. In other words, biased growth is a general feature of unmixed two-fluid simulations that manifests on scales where the individual transfer functions differ, consistent with the observations in Section 4.

We consider two simulations where CDM and baryon particles of uniform mass are evolved in an EdS universe ($\Omega_{c}=\Omega_{b}=0.5$; $h=1$): 1) a single-mode plane wave and 2) a three-mode perturbation. In both cases, the perturbations are aligned with the $x$-axis. Starting from two offset $128^{3}$ grids for baryons and CDM in an $L=32\Mpch$ box, we apply the ZA for an initial potential of the form $\phi^{\alpha}(q_{x})=\sum_{i}A^{\alpha}_{i}\cos(k_{i}q_{x})$, with $\alpha\in\{{c},{b}\}$ denoting CDM or baryons, and where $k_{i}=2\pi i/L$ for $i>0$ are the wave vectors. The amplitudes $A^{\alpha}_{i}$ are chosen so that the initial conditions at $z_{i}=200$ remain in the linear regime, i.e., well before shell-crossing occurs. For the single-mode simulation, we fix $A^{\alpha}_{i}=0$ except for $i=5$ and choose $A^{c}_{5}=6.1\times 10^{4}\;(\mathrm{km}/\mathrm{s})^{2}$. To imitate the effect of separate transfer functions, the baryon amplitude is set to $A^{b}=A^{c}/4$. For the second setup, we select the $i=(2,5,9)$ modes to be non-zero. We set $A^{c}$ to $(9.5,15,8.5)\times 10^{3}\;(\mathrm{km}/\mathrm{s})^{2}$ for the CDM perturbations and $A^{b}=(A^{c}_{2},A^{c}_{5}/4,A^{c}_{9})$ for the baryons, i.e., only the $i=5$ baryon mode has a different amplitude. We then compute the gravitational accelerations on each particle for $\Delta_{F}=\Delta_{P}$ and $\Delta_{F}=\Delta_{P}/4$ using the HACC-PM force solver in the first time step.

The particle configurations and acceleration measurements are illustrated in Fig. 12. Starting with the plane wave setup, we find that the particle accelerations for baryons are damped (i.e., closer to zero) on both sides of the density peak for the $\Delta_{F}=\Delta_{P}/4$ force resolution, whereas the CDM accelerations are boosted. This also manifests itself in Fourier space, where we find a decrease (increase) in the amplitude of the acceleration spectrum for baryons (CDM) at $k=5\times 2\pi/L$ when increasing the force resolution. Notably, for $\Delta_{F}=\Delta_{P}$, the CDM and baryon acceleration spectrum is equivalent, as is expected if both species describe continuous fluids. Discreteness errors appear when sampling the forces on a finer resolution, and in the case of this toy model, the opposite effect on CDM and baryons can be understood when looking at the relative position of CDM and baryon particles around the density peaks. Due to the different displacement amplitudes for the individual species, the CDM and baryon particles form “pairs”, where the CDM particle is always on the outside with respect to the density peak. If the force resolution is fine enough to resolve the 2-body interactions, this alignment will enhance (suppress) the inward acceleration for CDM (baryons), boosting (damping) the growth of that specific mode. This consistent arrangement occurs as long as the displacements are small, i.e., at high redshift. For larger perturbations, the particle ordering will eventually invert at some positions.

We show the same statistics for the three-mode configuration in the lower panel of Fig 12. We again find that the acceleration amplitudes for the individual species are coherently biased on the scale of the $k=5\times 2\pi/L$ mode when using the finer force resolution. This is the only mode that was initialized with different amplitudes for CDM and baryons. In the acceleration spectrum we also notice that the $k=9\times 2\pi/L$ mode is higher for both species. Since this is a species-independent feature, it is likely of a different origin than the asymmetric biased growth seen at $k=5\times 2\pi/L$. This was confirmed using a separate simulation where we used equal baryon and CDM amplitudes for all three modes, in which case we found the same difference in $P_{a_{x}}$ between force resolutions. Hence, the “bump” is merely caused by under-resolving this mode fully with $\Delta_{F}=\Delta_{P}$, and not a manifestation of discreteness errors in two-fluid simulations.

These toy models also provide us with some insight into how the different mitigation strategies (discussed in Section 5) circumvent the large-scale discreteness errors. Glass configurations such as used in B20 add noise to the pre-IC positions which destroys the consistent arrangement of baryon and CDM particle pairs. Using the total matter transfer function to calculate the particle displacements (H20) similarly prevents the formation of pairs by offsetting the two species coherently. If the simulation is initialized at lower redshift, the displacement amplitudes compared to the pre-IC interparticle separation is larger, allowing some CDM particles to “overtake” baryon particles and preventing a coherent acceleration bias on the scale of the $k$-mode where the initial amplitudes are different. Finally, softening or removing the short-range interspecies interactions removes the biased accelerations originating from two-body interactions.

We close this section by emphasizing that the two configurations presented are highly idealized, with specially selected modes and amplitudes used solely for illustrative purposes. Cosmological simulations are in general more complicated with biases introduced on all relevant modes, where interpretation of error growth is more difficult. That being said, the simple case-study does provide a consistent picture with the cosmological two-fluid simulation results, where discretization errors are confined to the scales at which the individual transfer functions are different.