# Novel Adaptive softening for collisionless N-body simulations: eliminating spurious haloes

Alexander Hobbs, Justin I. Read, Oscar Agertz, Francesca Iannuzzi, Chris Power

► To cite this version:

Alexander Hobbs, Justin I. Read, Oscar Agertz, Francesca Iannuzzi, Chris Power. Novel Adaptive softening for collisionless N-body simulations: eliminating spurious haloes. Monthly Notices of the Royal Astronomical Society, 2016, 458 (1), pp.468-479. 10.1093/mnras/stw251. hal-01440662

HAL Id: hal-01440662

https://hal.science/hal-01440662

Submitted on 9 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Novel Adaptive softening for collisionless N-body simulations: eliminating spurious haloes

Alexander Hobbs,<fnref n="1" />* Justin I. Read,<fnref n="2" /> Oscar Agertz,<fnref n="2" /> Francesca Iannuzzi<fnref n="3" /> and Chris Power<fnref n="4" />

Accepted 2016 January 28. Received 2016 January 28; in original form 2015 March 28

## ABSTRACT

We describe a Novel form of Adaptive softening (NOVA) for collisionless N-body simulations, implemented in the RAMSES adaptive mesh refinement code. In RAMSES – that we refer to as a ‘standard N-body method’ – cells are only split if they contain more than eight particles (a mass refinement criterion). Here, we introduce an additional criterion that the particle distribution within each cell be sufficiently isotropic, as measured by the ratio of the maximum to minimum eigenvalues of its moment of inertia tensor: $\eta = \lambda_{max}/\lambda_{min}$. In this way, collapse is only refined if it occurs along all three axes, ensuring that the softening $\epsilon$ is always of order twice the largest interparticle spacing in a cell. This more conservative force softening criterion is designed to minimize spurious two-body effects, while maintaining high force resolution in collapsed regions of the flow. We test NOVA using an antisymmetric perturbed plane wave collapse (‘Valinia’ test) before applying it to warm dark matter (WDM) simulations. For the Valinia test, we show that – unlike the standard N-body method – NOVA produces no numerical fragmentation while still being able to correctly capture fine caustics and shells around the collapsing regions. For the WDM simulations, we find that NOVA converges significantly more rapidly than standard N-body, producing little or no spurious haloes on small scales. We will use NOVA in forthcoming papers to study the issue of halo formation below the free-streaming mass $M_{fs}$; filament stability; and to obtain new constraints on the temperature of dark matter.

**Key words:** methods: numerical – dark matter – large-scale structure of Universe.

## 1 INTRODUCTION

The N-body method is widely used for modelling the non-linear growth of structure in the Universe (e.g. Dehnen & Read 2011; Kuhlen, Vogelsberger & Angulo 2012). For collisionless non-relativistic ‘cold’ dark matter (CDM), it has been shown to be remarkably accurate, producing robust results that are numerically well converged across different implementations (e.g. Heitmann et al. 2008; Springel et al. 2008; Stadel et al. 2009; Kim et al. 2014). Such simulations provide an excellent match to the observed large-scale structure in the Universe (e.g. Springel, Frenk & White 2006), though on smaller scales – where baryons likely play a role (e.g. Navarro, Eke & Frenk 1996; Read & Gilmore 2005; Mashchenko, Wadsley & Couchman 2008; Pontzen & Governato 2012) – there are known discrepancies (e.g. Flores & Primack 1994; Moore 1994; Klypin et al. 1999; Moore et al. 1999).

Despite the successes of the N-body method, since its inception there have been concerns about the effect of discreteness errors on numerical accuracy and convergence (e.g. Melott et al. 1997; Splinter et al. 1998; Binney 2004; Diemand et al. 2004; Wang & White 2007; Romeo et al. 2008; Joyce, Marcos & Baertschiger 2009). These arise because the dark matter fluid is represented by a set of ‘particles’, each with mass typically in the range $10^3–10^6 M_\odot$. To avoid spurious scattering between these particles, the force is softened, for example using ‘Plummer’ (Plummer 1915) softening:

$$F_{ij} = \frac{Gm^2(\mathbf{x}_j - \mathbf{x}_i)}{(\epsilon^2 + |\mathbf{x}_i - \mathbf{x}_j|^2)^{3/2}}, \quad (1)$$

where $F_{ij}$ is the force between two particles i and j at positions $\mathbf{x}_{i,j}$; G is Newton’s gravitational constant; and $\epsilon$ is the force softening. Equation (1) ensures that the force is clipped at a constant value as two particles approach one another. This reduces spurious two-body scattering, but it does not prevent numerical relaxation

<footnote n="1">¹Institute for Astronomy, ETH Zürich, CH-8093, Switzerland</footnote>

<footnote n="2">²Department of Physics, University of Surrey, Guildford, GU2 7XH Surrey, UK</footnote>

<footnote n="3">³Laboratoire d'astrophysique de Marseille, UMR 7326, F-13388, Marseille, France</footnote>

<footnote n="4">⁴The University of Western Australia, 35 Stirling Highway, CRAWLEY WA 6009, Australia</footnote>

<footnote>* E-mail: ahobbs@phys.ethz.ch</footnote>
from occurring; that can only be combated by raising the number of particles (e.g. Dehnen 2001; Power et al. 2003; Binney 2004; Diemand et al. 2004; Dehnen & Read 2011).

For CDM simulations, numerical convergence appears to be very good (e.g. Heitmann et al. 2008). However, discreteness errors may yet play a role when attempting to calculate power spectra, mass functions and higher order halo statistics at percent level accuracy, as will be required by next generation cosmological probes (e.g. Amendola et al. 2013; Reed et al. 2013). More problematic, however, are simulations that model a sharp cut-off in the initial power spectrum, as in *warm dark matter*<fnref n="1" /> (WDM; Avila-Reese et al. 2001; Bode et al. 2001), or exotic inflationary models (Zentner & Bullock 2003). The first WDM simulations appeared to find evidence of fragmentation – smaller haloes forming later than larger ones – as evidenced by a sharp upturn in the halo mass function (e.g. Avila-Reese et al. 2001; Bode et al. 2001). However, this has now been traced to the numerical fragmentation of filaments due to discreteness effects. This is particularly worrisome since the ‘spurious haloes’ that form via this process diminish with particle number only as $N^{1/3}$, leading to extremely slow convergence (Wang & White 2007).

The likely reason for the formation of ‘spurious haloes’ in WDM simulations was only recently elucidated by Hahn et al. (2013). Using a new method for evolving collisionless fluids – where they track the foliations of the 3D dark matter phase sheet – they find that the spurious haloes result from large anisotropic force errors. With a more accurate force (as calculated by their new method), the spurious haloes are much reduced, and the resulting filaments are smooth.

While Hahn et al. (2013) present an elegant alternative to N-body simulations, their method is numerically expensive. Since they are required to track the folding of the phase sheet, at the centres of dark matter haloes where there are many such foliations they formally require an ever-increasing number of simulation elements (Hahn & Angulo 2016); without such refinement, unphysical behaviour occurs in high-density regions. By contrast, a key strength of the ‘standard’ N-body method is that, since the equations of motion are derived from a Hamiltonian, the time-averaged expectation value of the energy of a particle will be correct even if its orbital phase is wrong<fnref n="2" /> (e.g. Dehnen & Read 2011). This likely explains the success of the N-body method even at rather modest N. For example, Sellwood (2008) find convergence for their disc simulations

already with $N \sim 10^5$, despite earlier calculations suggesting some $\sim 10^8$ particles would be required to properly resolve resonances (Weinberg & Katz 2007).

The above motivates considering whether the classic N-body method cannot be improved. Two recent works have attempted to ‘repair’ N-body simulations in post-processing by pruning spurious haloes. Schneider, Smith & Reed (2013) propose fitting a power law to the artificial haloes and subtracting them away, taking advantage of the fact that spurious haloes are more prevalent in overdense regions. By contrast, Lovell et al. (2014) suggest an algorithm where subhaloes are removed from the mass function if: (i) their ‘protohaloes’ are highly flattened; and/or (ii) the subhaloes are below a mass cut; and/or (iii) the subhaloes are not present in a higher resolution simulation of the same halo. In this paper, we consider instead a modified force softening criterion. This is designed to improve the anisotropic force errors that are at the root of the problem, leading to a more faithful N-body method in the first place.

This paper is organized as follows. In Section 2, we briefly review different strategies in the literature for force softening and we present our new Novel form of Adaptive softening – NOVA – designed to minimize spurious two-body effects. In Section 3, we describe our NOVA algorithm in detail and its implementation in the RAMSES code. In Section 4, we compare NOVA to standard RAMSES for an asymmetric plane wave test and 0.2 keV WDM simulations. We focus here on presenting the first results from NOVA for the density field; mass function; and dark matter halo density profiles. A detailed analysis of halo formation below the WDM ‘free-streaming’ mass (see Section 4.2.1); filament stability; and obtaining new constraints on the temperature of dark matter using NOVA will be presented in forthcoming publications. Finally, in Section 6 we present our conclusions.

## 2 FORCE SOFTENING

Since spurious haloes in WDM simulations appear to result from anisotropic force errors, this suggests that a good place to begin in improving the N-body method is to take a critical look at how the force softening $\epsilon$ is chosen.<fnref n="3" /> Ideally, $\epsilon$ should be as small as possible such that the maximum possible force resolution is obtained for a given numerical cost. However, too small and spurious forces will creep in, potentially spoiling numerical convergence. We consider three different force softening strategies here which are as follows.

(i) *Minimizing two-body effects (POWER)*: a popular rule-of-thumb, that has been carefully tested on CDM simulations, follows from ensuring that two body forces are small as compared to the mean field (Power et al. 2003):

$$ \frac{Gm^2}{\epsilon^2} \sim \frac{1}{\alpha^2} \frac{GmM}{r^2} \Rightarrow \epsilon \sim \alpha \frac{r}{\sqrt{N}}, \qquad (2) $$

where $\alpha = 4$ is an empirically derived parameter (Power et al. 2003); $N = M/m$ is the number of particles inside some characteristic radius $r$ (Power et al. 2003 use the virial radius $r_{200}$); and $m$ is the dark matter particle mass.

<footnote n="3">³ Note that it is equally important to select an appropriate timestep for the particles (e.g. Power et al. 2003; Zemp et al. 2007). However, this is true both in the standard N-body method and in the ‘folding phase sheet’ model of Hahn et al. (2013). This suggests that either the choice of $\epsilon$ is more crucial than the choice of timestep, or that the timestep criteria typically used in N-body simulations (e.g. Dehnen & Read 2011) are adequate.</footnote>

<footnote n="1">¹ In WDM, it is supposed that the dark matter is non-relativistic for a time after decoupling, leading to a suppression in the growth of structure on small scales and at early times (e.g. Avila-Reese et al. 2001; Bode, Ostriker & Turok 2001). Typically, this is modelled as an exponential cut-off in the initial power spectrum and indeed throughout this paper, where we refer to WDM simulations, this is what we mean. Fully self-consistent WDM models (for example, sterile neutrinos) have more complex model-dependent power spectra than this (e.g. Boyarsky, Ruchayskiy & Shaposhnikov 2009). Furthermore, for hot dark matter, it can also become important to model the primordial velocity dispersion of the dark matter particles. This has been attempted only a few times in the literature, most likely because of the computational cost involved. A proper treatment requires us to replace each ‘CDM’ particle in the initial conditions by $\sim 1000-10 000$ lighter particles in order to well sample the local velocity distribution function at each point in the flow (e.g. Avila-Reese et al. 2001; Hahn, Abel & Kaehler 2013). As far as the authors are aware, at the time of writing, such an expensive approach has never been attempted.</footnote>

<footnote n="2">² This is only strictly true if a symplectic time integrator is used. The Leapfrog integrator typically employed in cosmological simulations is symplectic, but only for fixed timesteps (e.g. Dehnen & Read 2011).</footnote>
(ii) *Minimizing force errors* (DEHNEN): an alternative approach is to minimize errors coming from biased forces that occur if $\epsilon$ is too large, and noise that occurs if $\epsilon$ is too small (Dehnen 2001). This leads to a well-defined optimal force softening that depends on the particular gravitational potential being simulated (and the choice of softening kernel). For small $\epsilon$ and large $N$, Dehnen (2001) derive an analytic estimate for Plummer softening that scales as

$$ \epsilon \propto N^{-0.73}; \quad \text{valid for } \epsilon \ll r; N \gtrsim 10^5. \qquad (3) $$

(iii) *Minimizing scatter between an ensemble of N-body realizations* (ROMEO): finally, following earlier work by Melott et al. (1997) and Splinter et al. (1998), Romeo et al. (2008) take a different approach. They run ensembles of the same cosmological N-body simulation, varying only the random number seed and the force softening. They argue that the scatter in results (as measured by various metrics like the halo mass function) has a term that is physical (cosmic variance), and a term that is numerical (discreteness noise). The optimum force softening should minimize the discreteness noise and therefore should minimize the scatter between different realizations of the same simulation. Using a novel wavelet analysis, they empirically derive

$$ \epsilon \sim 2d, \qquad (4) $$

where $d$ is the mean interparticle spacing.

While each of the above approaches to force softening is conceptually different, they all point to a rather similar conclusion: $\epsilon$ must be adaptive, varying both in space and time: $\epsilon \equiv \epsilon(\mathbf{x}, t)$. To see this, it is instructive to consider a simple toy model where the dark matter is spherically distributed with a power-law density profile:

$$ \rho \propto r^{-\gamma} \Rightarrow N(<r) \propto r^{(3-\gamma)}, \qquad (5) $$

where $N(<r)$ is the cumulative number of particles within $r$. (This equation is only strictly valid near the centre of dark matter haloes.) For this toy model, our three criteria give rather different results, but all point towards $\epsilon$ being some function of radius $r$ and therefore of the local density:

$$ \epsilon \propto \rho^{-\kappa}; \quad \kappa > 0, \qquad (6) $$

where

$$ \kappa = \begin{cases} (\gamma - 1)/2\gamma & \text{power} \\ 0.73(\gamma - 3)/\gamma & \text{dehnen} \\ 1/3 & \text{romeo.} \end{cases} $$

There are several interesting points to note from equation (6). First, notice that the POWER criterion actually amounts to fixed softening if $\gamma = 1$, as is the case for the centres of CDM haloes (Dubinski & Carlberg 1991). This may explain why fixed softening simulations have performed so surprisingly well despite the natural expectation that $\epsilon$ ought to be adaptive. If we are to adapt $\epsilon$, however, the POWER criterion becomes potentially pathological. For $\gamma < 1$, $\kappa < 0$ and in shallow dark matter cusps or cores, the softening would actually increase with density. The DEHNEN criterion fares better in this regard, being well behaved for all $\gamma < 3$; however, more work is required to generalize it to larger radii where the softening will be large and equation (3) is then no longer valid. For these reasons, we consider only the ROMEO criterion from here on.

The ROMEO softening has (perhaps inadvertently) been extensively explored in the literature. Mesh-based methods like the RAMSES code (Teyssier 2002) tie the softening to the local cell size which is naturally adaptive. Cells are split if they have greater than $N_c$ particles, typically chosen to be $N_c \sim 8$ to achieve, on average, one

<figure><img src="image_3.png" /><figcaption>**Figure 1.** A schematic view of the need for a modified adaptive force softening for N-body simulations. The simulation begins with the distribution locally very close to isotropic (left), with mean interparticle spacing $\Delta x \sim \Delta y$. However, as collapse proceeds first along the shortest axis (in this case the x-axis), we quickly move to a situation where $\Delta x \ll \Delta y$ (right). Standard adaptive softening schemes adapt purely on the local density. In this case, we would pick a softening $\epsilon \propto \Delta x \equiv \epsilon_x$, making the softening too small in the y-direction. This could lead to spurious clumping along the filament. In NOVA, we measure the local anisotropy and do not refine $\epsilon$ if the anisotropy is too high. In the example pictured, this would set our softening proportional to the longest local axis of the collapse – in this case $\epsilon_y \propto \Delta y$. The softening remains isotropic, but is more conservative than standard schemes in regions of high anisotropy.</figcaption></figure>

particle per cell after cell refinement. This amounts to a scheme where $\epsilon \propto 1/\rho_{\text{local}}^{1/3}$, exactly as in the ROMEO force softening. Similar schemes have also been explored in Tree N-body codes. There, since the equations of motion are derived from a Hamiltonian, it is possible to craft a density-adaptive $\epsilon$ method that is manifestly conservative (Price & Monaghan 2007). Iannuzzi & Dolag (2011) have recently implemented this in the GADGET code (Springel 2005), finding that it leads to results in excellent agreement with the fixed $\epsilon$ case, while giving greater resolution for similar numerical cost (see also Kawata et al. 2013). Their results suggest that with or without the conservative correction terms, the halo mass function converges; however, without the correction there is a substantial suppression of low-mass haloes mass as compared to both the conservation-corrected and fixed softening simulations. In Appendix A, we show that such conservative corrections do not solve ‘spurious halo’ problem in WDM simulations. In fact, since the correction terms appear as a purely attractive force that points along the density gradient, they make the spurious halo problem worse.

In this paper, we present a Novel form of Adaptive softening – NOVA – designed to minimize spurious two-body effects. Like the ROMEO softening, we tie the softening to the local interparticle spacing $\epsilon \sim 2d$. However, for the first time we account for the fact that in cosmological simulations, collapse is expected to be locally anisotropic (e.g. Zeldovich 1978). Since collapse proceeds most rapidly along the short axis of the flow, in the early stages of collapse there will always be three interparticle spacings aligned along the short c-, intermediate b-, and long a-axis (see Fig. 1). If we adapt $\epsilon$ purely on density, this amounts to an optimistic criterion $\epsilon \sim 2c$ that actually violates the ROMEO criterion along the long axis, leading to potentially large spurious two-body forces. Instead, our new NOVA method allows $\epsilon$ to be adapted on density only if the collapse is sufficiently isotropic. We show that this simple change to the N-body algorithm prevents the formation of ‘spurious’ haloes along filaments. [Note that, while in this paper we adapt on density, in principle NOVA can be applied to any adaptive softening scheme that obeys equation (6). We defer such generalizations to future work.]
# 3 NUMERICS

## 3.1 The RAMSES 'standard N-body' code

We carry out cosmological N-body simulations using the Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002). The collisionless dark matter dynamics are evolved using the particle-mesh technique (see e.g. Hockney & Eastwood 1988), with gravitational accelerations computed from the gravitational potential on the mesh. The gravitational potential is calculated by solving the Poisson equation using the multigrid method (Guillet & Teyssier 2011) for all refinement levels.

Note that while we describe RAMSES as a 'standard N-body' code, it actually differs from most Tree N-body codes in an important respect; the softening $\epsilon$ is automatically adapted according to the ROMEO criterion if the refinement strategy is based on reaching a critical number of particles per cell. Since this is true for all adaptive mesh refinement schemes in the literature to date, however, we still refer to this as 'standard'. We compare and contrast RAMSES with some GADGET Tree N-body simulations that use both fixed and adaptive softening in Appendix A.

## 3.2 The NOVA algorithm

In this section, we describe our Novel Adaptive force softening algorithm: NOVA. This is a modified cell splitting criterion implemented in the RAMSES code. Normally, cells are split if they contain more than some critical number of particles in a cell: $N_{cell} > N_c$. Here, we add an additional criterion that the cell is sufficiently isotropic as measured by its moment of inertia tensor

$$I_i = \sum_{j=1}^{N} \begin{bmatrix} (y_{ij}^2 + z_{ij}^2) m_j & -x_{ij} y_{ij} m_j & -x_{ij} z_{ij} m_j \\ -x_{ij} y_{ij} m_j & (x_{ij}^2 + z_{ij}^2) m_j & -y_{ij} z_{ij} m_j \\ -x_{ij} z_{ij} m_j & -y_{ij} z_{ij} m_j & (x_{ij}^2 + y_{ij}^2) m_j \end{bmatrix}, (7)$$

where $x_{ij} = x_i - x_j$, $y_{ij} = y_i - y_j$, $z_{ij} = z_i - z_j$.

We compute the eigenvalues of the matrix $I_i$: $\lambda_1, \lambda_2, \lambda_3$, which are sorted so that $\lambda_1 > \lambda_2 > \lambda_3$, and take the ratio $q_i \equiv \lambda_1/\lambda_3$ to be a measure of the (spatial) anisotropy in the particle distribution. We explored the effect of other measures of the local anisotropy: a simple cubic tensor,<fnref n="4" /> and the velocity dispersion tensor. Both also led to a reduction in the formation of spurious haloes as compared to the standard N-body method, but neither performed as well as our default scheme using the moment of inertia tensor.

Cells are split if $N_{cell} > N_c$ and $q_i < \eta$, where $\eta$ is a parameter that controls the amount of anisotropy allowed for splitting to occur. The effect of this is shown schematically in Fig. 1. By refining only where the particle distribution is locally isotropic, NOVA effectively picks the most conservative local force softening. As a result, the softening in anisotropic regions is always at least twice the longest interparticle spacing in a cell: $\epsilon \gtrsim 2 \max[d]$. This means that the softening is somewhat overestimated with respect to the short axis. However, too large softening only affects the computational efficiency (since for the same particle number we have reduced force resolution); whereas too small softening can – through two body effects – be much more problematic. An alternative way to think of the algorithm is that it does not refine unless collapse is occurring along all three axes. This typically reduces refinement in filamentary or elongated structures. Improved schemes could maintain the local

<table><caption>**Table 1.** Using cell particle noise to choose the anisotropy parameter $\eta$ for a given number of particles per cell $N_c$. The columns show $N_c$; the variance in anisotropy parameter for random realizations of a uniform particle distribution $\sigma_q$; and our choice of $\eta$ derived from this analysis: $\eta = 1 + 0.5\sigma_q$. Our approach has two desirable properties: (i) $\eta$ is set by the noise level for a cell; (ii) as a result, $\eta$ naturally shrinks with $N_c$. The grey row marks our default choice of parameters.</caption><thead><tr><th>N<sub>c</sub></th><th>σ<sub>q</sub></th><th>η</th></tr></thead><tbody><tr><td>8</td><td>0.748 626</td><td>1.374</td></tr><tr><td>16</td><td>0.304 557</td><td>1.152</td></tr><tr><td>32</td><td>0.169 258</td><td>1.084</td></tr><tr><td>64</td><td>0.102 809</td><td>1.051</td></tr><tr><td>128</td><td>0.066 1424</td><td>1.033</td></tr><tr><td>256</td><td>0.044 8105</td><td>1.022</td></tr><tr><td>512</td><td>0.030 1317</td><td>1.015</td></tr></tbody></table>

force resolution in filaments by allowing for anisotropic force softening as in earlier work by Shapiro et al. (1996). We will consider such schemes in a forthcoming paper.

Note that in practice each cell will contain many dark matter streams that are distinct in phase space, but overlapping in configuration space (e.g. Vogelsberger & White 2011). Recently, some authors (Hahn et al. 2013; Sousbie & Colombi 2015; Hahn & Angulo 2016) have attempted to actually resolve these streams, projecting them to determine the matter density and then solve the Poisson equation. NOVA similarly to the standard N-body method, makes the approximation that we can discard such stream information at approximately the force softening scale. Exploring the accuracy and convergence of structure formation simulations under this approximation is one of the goals of this work.

## 3.3 The choice of $\eta$ and $N_c$

In the limit $N \to \infty$, we would ideally have an anisotropy bound of $\eta = 1 - i.e. cell splitting is allowed only for purely isotropic cells. However, in practice, noise in the particle distribution makes it undesirable to set $\eta = 1$ exactly. Here, we choose as default $N_c = 32$ which is chosen to ensure that there are always enough particles in a cell that $I_i$ can be reliably measured ($N_c = 32$ ensures at least $\sim 3$ particles per spatial dimension); and we select $\eta = 1.08$. The latter number is chosen by drawing 32 particles at random from a uniform density distribution and calculating the distribution function of $q$. We choose $\eta = 1 + 0.5\sigma_q$, where $\sigma_q$ is the variance in $q$ for this random sampling; similar results for a range of $N_c$ are reported in Table 1. This has two desirable properties: (i) $\eta$ is set by the noise level for a cell; (ii) as a result, $\eta$ naturally shrinks with $N_c$. Note that the above implies that varying $\eta$ with the number of particles in a cell, or with the refinement level may give improved performance over the fixed $\eta$ scheme we explore here; such improvements are beyond the scope of this work. We explore the effect of varying $N_{cell}$ and $\eta$ in Section 4.3.

## 3.4 Numerical performance

As with any new numerical algorithm, we would be remiss not to discuss its numerical cost. We find that the relative performance of NOVA to RAMSES depends on the precise choice of parameters and

<footnote n="4">⁴ Of the form $X_{\alpha\beta} m_j$ where $1 \le \alpha \ge 3$, $1 \le \beta \ge 3$ refer to Cartesian basis vectors such that e.g. $X_{11} \equiv x_{ij}x_{ij}$.</footnote>
<table><caption>**Table 2.** The cosmological simulations and their parameters. The columns show from left to right: the simulation label; the particle resolution; the minimum number of particles per cell $N_c$; the anisotropy parameter $\eta$ (see Section 3); the dark matter particle mass $m_{\text{part}}$; and the maximum refinement level reached by that simulation.</caption><thead><tr><th>Label</th><th>Resolution</th><th>N<sub>c</sub></th><th>η</th><th>m<sub>part</sub> (10<sup>10</sup> M<sub>☉</sub> h<sup>-1</sup>)</th><th>Max. refinement level</th></tr></thead><tbody><tr><td>RAMSES-256</td><td>256<sup>3</sup></td><td>32</td><td>none</td><td>7.1 × 10<sup>-2</sup></td><td>14</td></tr><tr><td>RAMSES-512</td><td>512<sup>3</sup></td><td>32</td><td>none</td><td>8.9 × 10<sup>-3</sup></td><td>16</td></tr><tr><td>RAMSES-1024</td><td>1024<sup>3</sup></td><td>32</td><td>none</td><td>1.1 × 10<sup>-3</sup></td><td>17</td></tr><tr><td>NOVA-256</td><td>256<sup>3</sup></td><td>32</td><td>1.08</td><td>7.1 × 10<sup>-2</sup></td><td>13</td></tr><tr><td>NOVA-512</td><td>512<sup>3</sup></td><td>32</td><td>1.08</td><td>8.9 × 10<sup>-3</sup></td><td>15</td></tr><tr><td>NOVA-1024</td><td>1024<sup>3</sup></td><td>32</td><td>1.08</td><td>1.1 × 10<sup>-3</sup></td><td>16</td></tr><tr><td>NOVA-256-Nc128</td><td>256<sup>3</sup></td><td>128</td><td>1.03</td><td>7.1 × 10<sup>-2</sup></td><td>7</td></tr><tr><td>NOVA-512-Nc128</td><td>512<sup>3</sup></td><td>128</td><td>1.03</td><td>8.9 × 10<sup>-3</sup></td><td>8</td></tr><tr><td>NOVA-1024-Nc128</td><td>1024<sup>3</sup></td><td>128</td><td>1.03</td><td>1.1 × 10<sup>-3</sup></td><td>14</td></tr></tbody></table>

problem setup. Our default NOVA scheme is ≳10 times faster than RAMSES with $N_c = 8$. However, this is simply because NOVA refines less. When comparing with standard RAMSES with $N_c = 32$ (that refines only one level deeper than NOVA; see Table 2), NOVA is almost the same speed. If we compare NOVA and RAMSES for simulations where NOVA refines similarly to RAMSES, NOVA is about ~20 per cent slower. The fact that the numerical cost is primarily set by the maximum refinement level rather than the particle number means that when going to high resolution with lower refinement, it is memory usage that becomes a problem, rather than CPU.

# 4 RESULTS

## 4.1 The asymmetric plane wave (Valinia) test

We set up initial conditions for an asymmetric plane wave test as in Valinia et al. (1997) and Hahn et al. (2013). This simple 2D test allows us to make a check on two-body discreteness effects without running a full cosmological volume. The plane wave is setup along the x-direction with the following (sinusoidal) phase perturbation in the y-direction:

$$ \phi(\mathbf{x}) = \bar{\phi} \cos \left( k_p \left[ x + \epsilon_a \frac{k_p}{k_a^2} \cos k_a y \right] \right), \quad (8) $$

where $k_p = 2\pi/L$, $k_a = 4\pi/L$ and $\epsilon_a = 0.2$. L refers to the size of the simulation box. $\bar{\phi}$ sets the value of expansion factor at which the first shell crossing occurs – this is set to $a_c = 1/7.7 \simeq 0.13$. The initial particle positions and velocities were obtained by applying the Zel'dovich approximation (Zel'dovich 1970) to an unperturbed regular Cartesian lattice.

The plane wave is allowed to evolve under the pure gravitational potential of the particles up until $a = 1$. The results for our default choice of $N_c = 32$ and $\eta = 1.08$ are shown in Fig. 2. The left-hand panels show the particle distribution; the right the AMR refinement map. Notice that in RAMSES, the filaments break up into regularly spaced clumps. This occurs because the standard cell splitting criterion refines on the filaments (see top-right panel). By contrast, in NOVA no such refinement occurs (see bottom-right panel) and the filaments are smooth. Both algorithms refine on the bound structures that form at nodes, capturing the same caustic/shell-like structures reported in Hahn et al. (2013).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** $N = 512^3$ run for the asymmetrically perturbed plane wave (Valinia) test, comparing standard RAMSES (top) and our default NOVA method (bottom). The left-hand panels show the particle distribution; the right the AMR refinement map. Notice that in RAMSES, the filaments break up into regularly spaced clumps. This occurs because the standard cell-splitting criterion refines on the filaments (see top-right panel). By contrast, in NOVA no such refinement occurs (see bottom-right panel) and the filaments are smooth. Both algorithms refine on the bound structures that form at nodes, capturing the same caustic/shell-like structures reported in Hahn et al. (2013).</figcaption></figure>

## 4.2 WDM simulations

### 4.2.1 Initial conditions and simulation analysis

The WDM simulations were set up as in Power et al. (2003). We used cosmological parameters $\Omega_0 = 0.27$, $\Omega_\Lambda = 0.73$, $h = 0.705$ and $\sigma_8 = 0.81$ at $z = 0$ (Komatsu et al. 2011). Initial conditions were created by generating a statistical realization of a Gaussian random density field in Fourier space, with variance given by the linear matter power spectrum, and the Zel'dovich approximation used to compute initial particle positions and velocities. A CDM power spectrum was obtained by convolving the primordial power spectrum $P(k) \propto k^{n_{\text{spec}}}$ with the transfer function appropriate for our chosen set of cosmological parameters, computed using the Boltzmann code CAMB (see Lewis, Challinor & Lasenby 2000).
The power spectrum for the WDM model was then obtained a la Bobe et al. (2001), by filtering the CDM power spectrum with an additional transfer function of the form

$$T^{\text{WDM}}(k) = \left( \frac{P^{\text{WDM}}(k)}{P^{\text{CDM}}(k)} \right)^{1/2} = [1 + (\alpha k)^{2\nu}]^{-5/\nu}, \quad (9)$$

where $\alpha$ is a function of the WDM particle mass (equation 9 in Bode et al. 2001), $k$ is the wavenumber and $\nu = 1.2$ is a numerical constant.

We chose a WDM thermal relic mass of $m_\chi = 0.2$ keV. Following Schneider et al. (2012), this corresponds to an effective ‘free-streaming’ scale of

$$\lambda_{\text{fs}}^{\text{eff}} \approx 0.049 \left(\frac{m_{\chi}}{\text{keV}}\right)^{-1.11} \left(\frac{\Omega_{\chi}}{0.25}\right)^{0.11} \left(\frac{h}{0.7}\right)^{1.22} \text{Mpc h}^{-1} \quad (10)$$

which for $\Omega_\chi = 0.25$; $h = 0.73$; and $m_\chi = 0.2$ gives $\lambda_{\text{fs}}^{\text{eff}} = 0.308 \text{ Mpc h}^{-1}$. And a ‘free-streaming mass scale’

$$M_{\text{fs}} = \frac{4\pi}{3} \rho \left( \frac{\lambda_{\text{fs}}^{\text{eff}}}{2} \right)^3, \quad (11)$$

where $\rho$ is the mean background density of the Universe (that is a function of cosmology and redshift $z$; $\overline{\rho(z=0)} = 277.3h^2 M_{\odot} \text{ kpc}^{-3}$; Peacock 1999). For purely linear collapse with no mode-coupling or fragmentation, no haloes should form below $M_{\text{fs}}$. At redshift $z=0$, and assuming the above cosmological parameters and WDM thermal relic mass, we have $M_{\text{fs}} = 2.26 \times 10^9 M_{\odot}$.

A second scale of interest is the length-scale at which the WDM transfer function is reduced by half: the ‘half-mode’ length:

$$\lambda_{\text{hm}} \approx 13.93 \lambda_{\text{fs}}^{\text{eff}} \quad (12)$$

which also has an associated mass scale, the ‘half-mode mass’:

$$M_{\text{hm}} \approx 2.7 \times 10^3 \text{ M}_{\text{fs}}. \quad (13)$$

This is the mass scale at which we expect the WDM mass function to noticeably deviate from the CDM case.

It is not clear if haloes should exist below $M_{\text{fs}}$ in WDM structure formation. Angulo, Hahn & Abel (2013) use their new ‘folding phase sheet’ method to argue that there are no haloes below $M_{\text{fs}}$ (other than substructure haloes that originate from haloes more massive than $M_{\text{fs}}$). However, this result relies on some manual pruning of the halo mass function – required due to errors in the halo finding algorithm. The halo finder we employ is AHF (Gill, Knebe & Gibson 2004; Knollmann & Knebe 2009) which has been extensively tested in CDM simulations and shown to be among the most robust in the field. However, we verify that halo finding in WDM is indeed a thorny issue, and discuss NOVA results for haloes below $M_{\text{fs}}$ in Section 4.2.2.

Note that we deliberately choose a small $m_\chi = 0.2$ keV even though such a low thermal relic mass is already ruled out by constraints from the Lyman-$\alpha$ forest (e.g. Viel et al. 2013). The reason for this is that it ensures that WDM effects will appear at large mass, making it computationally efficient to test our methodology (cf. Hahn et al. 2013). We will present NOVA simulations of particle masses close to the current observational constraints (and comparisons with data) in future work. A full list of all simulations run and their parameters is given in Table 2.

### 4.2.2 The density field

Fig. 3 shows a full box view of the projected dark matter density field in RAMSES and NOVA. Bound haloes are identified using AHF

(see Section 4.2.1); these are overplotted as red filled circles, with a size proportional to their virial radii. In the RAMSES simulation (left-hand panels), the filaments break up into many small haloes. These ‘spurious’ haloes have already been shown to be numerical artefacts (e.g. Wang & White 2007) and it is encouraging that they are gone in NOVA (bottom panels). Once these structures are removed, NOVA does a good job of capturing the caustics, fine shells, and criss-crossing filaments that surround galaxies.

Where caustic and shell structures overlap, AHF often identifies ‘haloes’, yet it is not clear if such structures are real or simply transient. We call these ‘caustic’ haloes. Their existence – if real – is important as it implies that haloes can indeed form below the free-streaming mass $M_{\text{fs}}$. A detailed exploration of this requires improving the halo finder and studying these structures carefully as a function of time. We will explore this in detail in a forthcoming publication.

### 4.2.3 The halo mass function and convergence

In Fig. 4, we compare the AHF cumulative halo mass function for RAMSES (blue) with NOVA (red) at three numerical resolutions: $N = 256^3$ (dashed), $512^3$ (dotted) and $1024^3$ (solid) particles. Overplotted is the cumulative mass function for the equivalent CDM simulation (black), and the free streaming mass $M_{\text{fs}}$ at redshift $z=0$ (dotted line).

The RAMSES simulations show very poor numerical convergence, with a prominent upturn that shifts extremely slowly to lower mass as the resolution is increased. It is interesting that our results for mass function convergence in RAMSES are somewhat worse than reported in Wang & White (2007). This likely owes to the fact that we use a different halo finder; we find that switching off the ‘unbind’ feature in AHF leads to very different behaviour, illustrating how sensitive results for WDM simulations are to the choice of halo finder and its chosen settings.

By contrast, NOVA shows much better numerical convergence. The mass function rises slowly at the low-mass end with resolution, while remaining unchanged at high mass. Below the free streaming scale, we find nearly an order of magnitude fewer haloes in NOVA than in RAMSES. We do however find a tendency for the lowest resolutions to oversuppress haloes at intermediate mass, although this goes away quickly with increasing resolution.

Our chosen halo finder AHF is known to perform well on CDM simulations (Knebe et al. 2011), but has not been used on ‘spurious halo free’ WDM simulations before. As discussed in Section 4.2.2 (and see Fig. 3), for the simulations we present here, AHF almost certainly misidentifies some features as haloes. For this reason, we defer a more careful analysis of haloes below the free streaming mass $M_{\text{fs}}$ to future work. We discuss how our results compare with other recent determinations in the literature in Section 5.

Note that our NOVA algorithm does lead to increasing refinement with resolution. This is shown in the final column of Table 2, where we list the maximum refinement level reached for all WDM simulations in this paper. The NOVA simulations typically reach ~1 refinement level less than standard RAMSES, but none the less they do continue to refine with increasing resolution. We discuss the effect of varying $\eta$ on the maximum refinement reached and on numerical convergence in Sections 4.3 and 5.1.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** A comparison of *RAMSES* and *NOVA* for a 0.2 keV WDM simulation with $N = 512^3$ particles. The top panels show the full 50 Mpc $h^{-1}$ box; the bottom panels highlight a zoomed in region, as marked. Haloes identified using the halo finder (AHF) are marked by the red filled circles; their size is proportional to their virial radii. Notice that in *RAMSES*, many small and regularly spaced haloes – ‘spurious haloes’ – form along the filaments; in *NOVA* these are gone. Notice further that in the zoom panel for *NOVA*, there are several cases of haloes identified by AHF that may not correspond to genuine bound structures that could host galaxies. We highlight two of these ‘caustic’ haloes as examples. Some of these haloes likely owe to overlapping caustics in the WDM density field and would not host galaxies; others, however, may be genuine non-linear structures that form at an overlap between caustics or as filaments intersect. This latter possibility is very interesting as it would imply that haloes can form below the free streaming mass $M_{fs}$ in WDM. We will explore this further in a forthcoming paper.</figcaption></figure>

### 4.2.4 Halo density profiles

In Fig. 5, we show dark matter halo density profiles in the highest resolution (1024³) *RAMSES* and all three resolutions in our fiducial $N_c = 32$ *NOVA* simulations (left-hand panel) for both low-mass (~10¹² M☉) and high-mass (~10¹⁴ M☉) haloes. It has already been reported previously in the literature that simply forbidding

refinement will reduce spurious two-body effects in collisionless N-body simulations (e.g. Melott et al. 1997; Splinter et al. 1998; Hahn et al. 2013). However, reducing refinement everywhere also leads to the centres of haloes – where galaxies actually reside – becoming unresolved. In *NOVA*, we attempt to obtain the best of both worlds as the algorithm leads to derefinement in highly anisotropic
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dark matter cumulative halo mass functions in RAMSES (blue) and NOVA (red) for a 0.2 keV WDM simulation. In both cases, three different resolutions are marked: $N = 256^3$ (dashed); $N = 512^3$ (dotted); and $N = 1024^3$ (solid) particles. Overlaid is the curve expected for the equivalent CDM simulation (black). Notice also that the standard RAMSES simulations converge very slowly with increasing resolution, showing a characteristic upturn at low mass that only strengthens as the resolution increases. By contrast, NOVA converges rapidly 'from below'. At $N = 1024^3$ particles, the number of haloes at low mass in NOVA is suppressed with respect to standard RAMSES by nearly an order of magnitude. Finally, notice that even in NOVA, the cumulative mass function does not reach a plateau below the free-streaming mass scale: there are significant numbers of haloes below $M_{\text{fs}}$ (vertical dotted line; see Section 4.2.1). See the text for further discussion of this.</figcaption></figure>

regions while still refining on halo centres (see Fig. 2, bottom-right panel). For this reason, the NOVA haloes at our highest resolution have a density profile in excellent agreement with the equivalent RAMSES simulation. The ~$10^{12}$ M$_\odot$ mass haloes in NOVA are slightly shallower in the two lower resolutions ($256^3$, $512^3$) than their RAMSES counterparts, reflecting the lower refinement level reached. The effect diminishes with resolution, however, demonstrating that NOVA is convergent. In the right-hand panel of Fig. 5, we plot the largest ($10^{14}$ M$_\odot$) for alternative NOVA runs with $N_c = 128$, again comparing with the highest resolution fiducial $N_c = 32$ RAMSES simulation. We see a more extreme effect of the lower refinement, where the lowest two resolutions are substantially underresolved in the halo centre due to a ~7 levels difference in refinement compared to the $1024^3$ run.

## 4.3 Varying $N_c$ and $\eta$: should some filaments fragment after all?

In this section, we study the effect of varying the minimum number of particles in a cell $N_c$ and the anisotropy parameter $\eta$. For this paper, we use throughout $\eta = 1 + 0.5\sigma_\eta$ as outlined in Section 3. We defer a detailed analysis of $\eta$ – in particular allowing a time or spatially varying $\eta$ – to future work.

In Fig. 6, we show a zoom-in on a particularly massive filament that can be seen in the top middle of the full simulation box shown in Fig. 3. We focus on this particular filament because we find that it is very hard to avoid it fragmenting. At low resolution ($256^3$; left-hand panels), the filament is completely smooth. However, for our default choice of $N_c = 32$, already at $N = 512^3$, the filament begins to fragment (top middle panel), as does another massive filament that connects to the large halo at the edge of the box (see yellow circles). As we raise the numerical resolution $N$ at fixed $N_c$ and $\eta$, the filament fragments further (top-right panel). If we raise $N_c$ and lower $\eta$ according to our ‘noise criteria’ (Section 3; bottom panels) then at $N = 512^3$ the filament becomes once again smooth (bottom middle panel). However, even for these NOVA parameters, raising the resolution to $N = 1024^3$ results in the filament fragmenting (bottom-right panel). Indeed, the filament breaks up into structure

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left-hand panel: dark matter density profiles in a 0.2 keV WDM simulation for example haloes with virial masses ~$10^{14}$ M$_\odot$ (upper lines) and ~$10^{12}$ M$_\odot$ (lower lines) in RAMSES (blue) and NOVA (red) in the $256^3$ (dotted), $512^3$ (dashed) and $1024^3$ (solid) resolution simulations. The profiles correspond to the same haloes at each resolution in each method. Right-hand panel: the example ~$10^{14}$ M$_\odot$ halo in the $N_c = 128$ NOVA runs (all three resolutions; red) compared to the same halo in the $N_c = 32$ RAMSES run (blue). The large discrepancy between the ($256^3$, $512^3$) compared to the $1024^3$ in 128 per cell NOVA is due to substantial differences in the maximum refinement level: in this NOVA run the two lowest resolutions refine only down to levels (7, 8), respectively, while the highest resolution goes all the way to level 14, a factor of (128, 64) increase in effective resolution over the other two in the refined region (the centre of the halo).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Zoom-in on a particularly massive filament (taken from the top middle of the full simulation box; see Fig. 3) for varying $N_c$; $\eta$; and numerical resolution $N$, as marked on each panel. We focus on this particular filament because we find that it is very hard to avoid it fragmenting. At low resolution ($256^3$; left-hand panels), the filament is completely smooth. However, for our default choice of $N_c = 32$, already at $N = 512^3$, the filament begins to fragment (top middle panel), as does another massive filament that connects to the large halo at the edge of the box (see yellow circles). As we raise the numerical resolution $N$ at fixed $N_c$ and $\eta$, the filament fragments further (top-right panel). If we raise $N_c$ and lower $\eta$ according to our ‘noise criteria’ (Section 3; bottom panels) then at $N = 512^3$ the filament becomes once again smooth (bottom middle panel). However, even for these NOVA parameters, raising the resolution to $N = 1024^3$ results in the filament fragmenting (bottom-right panel). Interestingly, the largest structures that form in filaments appear in the same place regardless of our choice of NOVA parameters; these are marked by the yellow circles.</figcaption></figure>

in any version of NOVA where we permit refinement within the filament. Interestingly, the largest structures that form both in this filament and the one to its top left appear always in the same place regardless of our choice of NOVA parameters. These structures are marked by the yellow circles (yellow is chosen to avoid confusion with AHF haloes in previous plots that are marked in red). The fact that such structures are challenging to avoid and yet also appear always at the same locations suggests that they may be physically correct. Such behaviour is certainly rather different from the regularly spaced fragments that form along filaments in the standard RAMSES simulations (Fig. 3). We discuss this further in Section 5.

# 5 DISCUSSION

## 5.1 Numerical convergence in NOVA

It is challenging to quantitatively test convergence in NOVA because of uncertainties in the halo finding. For example, Wang & White (2007) using Friends-of-Friends (FoF) report a definite, if slow, shift of the upturn in the halo mass function to lower mass with increasing resolution. Using AHF, we find also a shift in the upturn but it is substantially slower (see Fig. 4). Taken at face value, it is

not at all clear that RAMSES will converge on the correct solution with increasing $N$.

By contrast, convergence in NOVA seems healthier. At low mass, the mass function is suppressed causing a ‘convergence from below’ with increasing $N$. At low resolution, there is a clear oversuppression in haloes even at high mass (at $N = 256^3$ RAMSES outperforms NOVA). But with increasing $N$, the situation rapidly improves. We find a gently rising mass function with few low-mass haloes, but none the less some. Again, owing to difficulties with halo finding, we cannot yet determine whether these low-mass haloes are real; this is either a shortcoming of the halo finder, or evidence that NOVA requires further improvement. We defer a careful analysis of this to a future work where we will improve on the halo finder and study the time evolution of haloes.

## 5.2 Comparison with recent work in the literature

Despite the nearly an order of magnitude suppression in spurious haloes in NOVA, at first sight the mass function appears to rise substantially more steeply towards low mass than that reported recently in Angulo et al. (2013); using their new ‘folding phase sheet’ methodology, they find no haloes below $M_{\text{fs}}$. However, there are
three confounding factors that make a comparison difficult. First, Angulo et al. (2013) use a FoF halo finder similarly to Wang & White (2007), whereas we use AHF. Secondly, they throw out all haloes with overlapping virial radii. This has the effect that all sub-structure haloes are removed, providing a better comparison with semi-analytic theories for structure formation like Press–Schechter (Press & Schechter 1974). However, such a method will also remove half of all binary systems. Finally, they further prune their mass function through a visual inspection of haloes. This is necessary due to the FoF algorithm picking up many false positives. We find fewer false positives when using AHF, but many of the issues they report with halo finding in WDM resonate with our findings here. As pointed out in Section 4.2.2 (and see Fig. 3), we do see many suspicious structures identified by AHF as bound haloes that are unlikely to host galaxies. Given the difficulty of halo finding in WDM, we defer a more careful analysis of haloes below $M_{\text{fs}}$ to a future work where we will look critically at the time evolution of haloes in NOVA and attempt to improve on existing halo finding algorithms that have all been tuned to work well only on CDM simulations (e.g. Knebe et al. 2011).

Finally, given that the spurious structure appears to be a direct result of choosing the wrong softening length in highly anisotropic regions, it is tempting to suggest a form of anisotropic kernel as the solution (e.g. Shapiro et al. 1996; Romeo 1997, 1998). Such schemes come with their own complications, however. We will consider these in a forthcoming paper.

## 5.3 Filaments that physically fragment?

In crafting a scheme that avoids spurious fragmentation, we should be mindful of not throwing the baby out with the proverbial bathwater. Are we really sure that no filament should fragment, or that all such fragments reported in N-body simulations are spurious? Our results in Section 4.3 strongly suggest that at least some filaments may be unstable to physical fragmentation. From a theoretical point of view, this should perhaps not be surprising. It has long been known that infinite self-gravitating cylinders are unstable to fragmentation (Ostriker 1964; Fridman & Poliachenko 1984). However, cosmological filaments are both finite and expanding, and given sufficient expansion, self-gravitating cylinders become unconditionally stable (Schneider & Moore 2011). We will explore this further in future work, but the issue is an important one: if filaments are physically unstable then haloes can collapse below $M_{\text{fs}}$ and isolated small galaxies can be expected to exist even in warm/hot dark matter cosmologies.

## 6 CONCLUSIONS

We have introduced a Novel form of Adaptive softening (NOVA) for collisionless N-body simulations, implemented in the RAMSES adaptive mesh refinement code. In RAMSES – that we refer to as a ‘standard N-body method’ – cells are only split if they contain more than eight particles (a mass refinement criterion). We introduced an additional criterion that the cell be sufficiently isotropic, as measured by the ratio of the maximum to minimum eigenvalues of its moment of inertia tensor: $\eta = \lambda_{\text{max}}/\lambda_{\text{min}}$. In this way, collapse is only refined if it occurs along all three axes, ensuring that the softening $\epsilon$ is always of order twice the largest interparticle spacing in a cell. This more conservative force softening criterion was designed to minimize spurious two-body effects, while maintaining high force resolution in collapsed regions of the flow.

We tested NOVA using an antisymmetric perturbed plane wave collapse (‘Valinia’ test) before applying it to WDM simulations. Our key results are as follows.

(i) We used the Valinia test to show that – unlike the standard N-body method (RAMSES) – NOVA produces no numerical fragmentation while still being able to correctly capture high-density features like the fine caustics and shells around the collapsing regions.

(ii) For the WDM simulations, we found that NOVA converges significantly more rapidly than RAMSES, producing little or no spurious haloes on small scales. NOVA produces nearly an order of magnitude fewer dark matter haloes at low mass as compared to RAMSES, while still being able to correctly resolve high-density regions at the centres of massive haloes.

(iii) Despite the large reduction in low-mass haloes in NOVA, we still found haloes below the free streaming mass scale $M_{\text{fs}}$. Furthermore, these haloes increase in number (albeit slowly) as we increase the numerical resolution. Some of these likely owe to our halo finder (AHF) incorrectly labelling caustics and criss-crossing filaments as haloes. Others form as larger haloes that form above $M_{\text{fs}}$ are tidally stripped. However, some isolated low-mass structures appear to be real. Due to the difficulty of accurate halo identification in WDM, we defer a quantitative analysis of haloes below $M_{\text{fs}}$ to a forthcoming publication.

(iv) We highlighted two particularly massive filaments that fragment in any version of NOVA where refinement is allowed. Since the most massive fragments appear always at the same locations, we argue that these may be physical. We noted that infinite self-gravitating cylinders are unstable to collapse and so particularly massive cosmological filaments may be physically unstable. We will explore this further in future work, but the issue is an important one: if filaments are physically unstable then haloes can collapse below $M_{\text{fs}}$ and isolated small galaxies can be expected to exist even in warm/hot dark matter cosmologies.

We will use NOVA in forthcoming papers to study the issue of halo formation below $M_{\text{fs}}$; filament stability; and to obtain new constraints on the temperature of dark matter.

## ACKNOWLEDGEMENTS

JIR would like to acknowledge support from SNF grant PP00P2_128540/1. We would like to thank Oliver Hahn, Aurel Schneider, Romain Teyssier and Joachim Stadel for useful discussions. This research used the MÖNCH supercomputer owned by ETH Zürich, and the SURREY GALAXY FACTORY machine at the University of Surrey.

## REFERENCES

* Amendola L. et al., 2013, Liv. Rev. Relativ., 16, 6

* Angulo R. E., Hahn O., Abel T., 2013, MNRAS, 434, 3337

* Avila-Reese V., Colín P., Valenzuela O., D'Onghia E., Firmani C., 2001, ApJ, 559, 516

* Binney J., 2004, MNRAS, 350, 939

* Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93

* Boyarsky A., Ruchayskiy O., Shaposhnikov M., 2009, Ann. Rev. Nucl. Part. Sci., 59, 191

* Dehnen W., 2001, MNRAS, 324, 273

* Dehnen W., Read J. I., 2011, Eur. Phys. J. Plus, 126, 55

* Diemand J., Moore B., Stadel J., Kazantzidis S., 2004, MNRAS, 348, 977

* Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496

* Flores R. A., Primack J. R., 1994, ApJ, 427, L1
Fridman A. M., Poliachenko V. L., 1984, Physics of Gravitating Systems. II- Nonlinear Collective Processes: Non-linear Waves, Solitons, Collision- less Shocks, Turbulence, Astrophysical Applications. Springer-Verlag, New York

Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

Hahn O., Angulo R. E., 2016, MNRAS, 455, 1115

Hahn O., Abel T., Kaehler R., 2013, MNRAS, 434, 1171

Heitmann K. et al., 2008, Comput. Sci. Discovery, 1, 015003

Hockney R. W., Eastwood J. W., 1988, Computer Simulation using Particles. Hilger, Bristol

Iannuzzi F., Dolag K., 2011, MNRAS, 417, 2846

Joyce M., Marcos B., Baertschiger T., 2009, MNRAS, 394, 751

Kawata D., Okamoto T., Gibson B. K., Barnes D. J., Cen R., 2013, MNRAS, 428, 1968

Kim J.-h. et al., 2014, ApJS, 210, 14

Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82

Knebe A. et al., 2011, MNRAS, 415, 2293

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Komatsu E. et al., 2011, ApJS, 192, 18

Kuhlen M., Vogelsberger M., Angulo R., 2012, Phys. Dark Universe, 1, 50

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014, MNRAS, 439, 300

Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174

Melott A. L., Shandarin S. F., Splinter R. J., Suto Y., 1997, ApJ, 479, L79

Moore B., 1994, Nature, 370, 629

Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19

Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72

Ostriker J., 1964, ApJ, 140, 1529

Peacock J. A., 1999, Cosmological Physics. Cambridge Univ. Press, Cam- bridge

Plummer H. C., 1915, MNRAS, 76, 107

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Press W. H., Schechter P., 1974, ApJ, 187, 425

Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Reed D. S., Smith R. E., Potter D., Schneider A., Stadel J., Moore B., 2013, MNRAS, 431, 1866

Romeo A. B., 1997, A&A, 324, 523

Romeo A. B., 1998, A&A, 335, 922

Romeo A. B., Agertz O., Moore B., Stadel J., 2008, ApJ, 686, 1

Schneider A., Moore B., 2011, MNRAS, 415, 1569

Schneider A., Smith R. E., Macciò A. V., Moore B., 2012, MNRAS, 424, 684

Schneider A., Smith R. E., Reed D., 2013, MNRAS, 433, 1573

Sellwood J. A., 2008, ApJ, 679, 379

Shapiro P. R., Martel H., Villumsen J. V., Owen J. M., 1996, ApJS, 103, 269

Sousbie T., Colombi S., 2015, preprint (arXiv:1509.07720)

Splinter R. J., Melott A. L., Shandarin S. F., Suto Y., 1998, ApJ, 497, 38

Springel V., 2005, MNRAS, 364, 1105

Springel V., Frenk C. S., White S. D. M., 2006, Nature, 440, 1137

Springel V. et al., 2008, MNRAS, 391, 1685

Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, 398, L21

Teyssier R., 2002, A&A, 385, 337

Valinia A., Shapiro P. R., Martel H., Vishniac E. T., 1997, ApJ, 479, 46

Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D, 88, 4, 043502

Vogelsberger M., White S. D. M., 2011, MNRAS, 413, 1419

Wang J., White S. D. M., 2007, MNRAS, 380, 93

Weinberg M. D., Katz N., 2007, MNRAS, 375, 425

Zel'dovich Y. B., 1970, A&A, 5, 84

Zeldovich I. B., 1978, in Longair M. S., Einasto J., eds, Proc. IAU Symp. 79, Large-Scale Structures in the Universe. Reidel, Dordrecht, p. 409

Zemp M., Stadel J., Moore B., Carollo C. M., 2007, MNRAS, 376, 273

Zentner A. R., Bullock J. S., 2003, ApJ, 598, 49

APPENDIX A: WHY CONSERVATIVE ADAPTIVE SOFTENING DOES NOT SOLVE THE PROBLEM OF SPURIOUS HALOES

One of our earlier ideas for solving the spurious halo problem was to use conservative adaptive force softening as originally suggested by Price & Monaghan (2007) and implemented recently in the GADGET code (Springel 2005) by Iannuzzi & Dolag (2011). In this appendix, for completeness, we report the results of this experiment and explain why it failed to produce the desired effect. In fact, such conservative corrections to adaptive force softening make the problem of spurious haloes even worse!

As detailed in Price & Monaghan (2007), if the force softening varies in space and time then we can still construct a fully conservative N-body method by deriving the equations of motion from a discretized Hamiltonian. This results in an additional corrective force to the usual N-body equation of motion:

$$ \frac{d\mathbf{v}_{i,c}}{dt} = -\frac{G}{2} \sum_j m_j \left[ \frac{\xi_i}{\Omega_i} \frac{\partial W_{ij}(h_i)}{\partial \mathbf{r}_i} + \frac{\xi_j}{\Omega_j} \frac{\partial W_{ij}(h_j)}{\partial \mathbf{r}_j} \right], \quad (A1) $$

where

$$ \Omega_i \equiv \left[ 1 - \frac{\partial h_i}{\partial \rho_i} \sum_j m_j \frac{\partial W_{ij}(h_i)}{\partial h_i} \right] \quad (A2) $$

<figure><img src="image_9.png" /><figcaption>Figure A1. Cumulative mass function for GADGET adaptive softening runs, with (magenta) and without (brown) the conservative correction term, compared to the fixed softening reference run (black). Notice that the upturn is more pronounced when the correction term is used. These runs used $N = 128^3$ particles.</figcaption></figure>
and

$$ \xi_i \equiv \frac{\partial h_i}{\partial \rho_i} \sum_j \frac{\partial \phi_{ij}(h_i)}{\partial h_i}, \qquad (A3) $$

where G is Newton's gravitational constant; $m_j$ is the mass of particle j; $W_{ij}$ is a positive definite spherical smoothing kernel; $h_i$ is the smoothing length of particle j, here equated with the softening $h_j = \epsilon_j$; $\mathbf{r}_j$ is the position of particle j; and $\phi_{ij} \equiv \phi(|\mathbf{r}_i - \mathbf{r}_j|)$ is related to the gravitational potential between particle pairs (Price & Monaghan 2007).

The key thing to note from equation (A1) is that, since the kernel is positive definite, $\frac{\partial W_{ij}(h_i)}{\partial h_i}$ and similar terms are negative definite; the correction terms $\Omega_i$ are always of order unity; while the $\xi_i$ terms are also negative definite. This means that the force correction in equation (A1) is negative definite and will lead always to an increased gravitational force. This force will point to leading order along the density gradient:

$$ \nabla_i \rho_i = \sum_j m_j \nabla_i W_{ij}(h_i). \qquad (A4) $$

This is the trouble with the conservative correction terms. When the flow becomes highly anisotropic, the correction terms will act to increase the force along the short axis leading to even more artificial clumping (see Fig. 1).

In Fig. A1, we show the cumulative mass function for a 50 Mpc $h^{-1}$ 0.2 keV WDM simulation run using the standard GADGET N-body code (black); with adaptive force softening (brown); and with conservative adaptive force softening (magenta). For this comparison, we used a kernel neighbour number of $N_{neigh} = 60$ since this was shown to be optimal in Iannuzzi & Dolag (2011). Notice that the sharp upturn in the mass function (present in all simulations) is more prominent in the simulation with the conservative correction term (magenta) than even in the fixed softening case (black).

This paper has been typeset from a TeX/Latex file prepared by the author.