Numerical convergence of simulations of galaxy formation: the abundance and internal structure of cold dark matter haloes

TL;DR

This paper investigates how numerical parameters affect the properties of cold dark matter haloes in cosmological simulations, providing criteria for convergence and assessing the robustness of simulated halo structures.

Contribution

It offers new analytic estimates for the convergence radius in dark matter halo simulations based on particle number and softening, improving reliability of simulation results.

Findings

Median circular velocity profiles converge within 10% beyond the convergence radius.

Convergence depends on the relaxation timescale and particle number, with better results at higher values.

Analytic criteria relate the convergence radius to inter-particle spacing, aiding in simulation assessment.

Abstract

We study the impact of numerical parameters on the properties of cold dark matter haloes formed in collisionless cosmological simulations. We quantify convergence in the median spherically-averaged circular velocity profiles for haloes of widely varying particle number, as well as in the statistics of their structural scaling relations and mass functions. In agreement with prior work focused on single haloes, our results suggest that cosmological simulations yield robust halo properties for a wide range of gravitational softening parameters, $\epsilon$, provided: 1) $\epsilon$ is not larger than a "convergence radius", $r_{\rm conv}$, which is dictated by 2-body relaxation and determined by particle number, and 2) a sufficient number of timesteps are taken to accurately resolve particle orbits with short dynamical times. Provided these conditions are met, median circular velocity profiles converge to within $\approx 10$ per cent for radii beyond which the local 2-body relaxation timescale exceeds the Hubble time by a factor $\kappa\equiv t_{\rm relax}/t_{\rm H}\gt 0.177$, with better convergence attained for higher $\kappa$. We provide analytic estimates of $r_{\rm conv}$ that build on previous attempts in two ways: first, by highlighting its explicit (but weak) softening-dependence and, second, by providing a simpler criterion in which $r_{\rm conv}$ is determined entirely by the mean inter-particle spacing, $l$; for example, better than $10$ per cent convergence in circular velocity for $r\gt 0.05\,l$. We show how these analytic criteria can be used to assess convergence in structural scaling relations for dark matter haloes as a function of their mass or maximum circular speed.

Figures

Tables

Table 1. Table 1: Numerical aspects of our dark matter only simulations. The first column provides a run label, adopting the same nomenclature as Schaye et al ( 2015 ) . N p subscript 𝑁 p N_{\rm p} is the total number of simulation particles of mass m DM subscript 𝑚 DM m_{\rm DM} ; ϵ cm subscript italic-ϵ cm \epsilon_{\rm cm} and ϵ phys subscript italic-ϵ phys \epsilon_{\rm phys} the co-moving and maximum physical softening lengths, respectively; ϵ i / l subscript italic-ϵ 𝑖 𝑙 \epsilon_{i}/l are the softening lengths expressed in units of the mean inter-particle separation, l = L b / N p 1 / 3 𝑙 subscript L b superscript subscript 𝑁 p 1 3 l={\rm L_{b}}/N_{\rm p}^{1/3} ; f ϵ subscript 𝑓 italic-ϵ f_{\epsilon} is the softening length in units of the “fiducial” eagle value for a given m DM subscript 𝑚 DM m_{\rm DM} ; ErrTolIntAcc is gadget ’s integration accuracy parameter; z phys subscript 𝑧 phys z_{\rm phys} the redshift below which the softening remains fixed in physical (rather than comoving) units.

	Name	$m_{\mathrm{DM}}$	$N_p$	${ϵ}_{\mathrm{phys}}$	${ϵ}_{\mathrm{phys}}/l$	${ϵ}_{\mathrm{cm}}$	${ϵ}_{\mathrm{cm}}/l$	$f_{ϵ}$	ErrTolIntAcc	$z_{\mathrm{phys}}$
		[${10}^5{M}_{\odot }$]		$\left[\mathrm{pc}\right]$	$\left[{10}^{-2}\right]$	$\left[\mathrm{pc}\right]$	$\left[{10}^{-2}\right]$	$\left[ϵ/{ϵ}_{\mathrm{fid}}\right]$
	L0012N0752	1.8	7523	175.0	1.05	665.0	4.00	1.0000	0.025	2.8
	L0012N0752	1.8	7523	43.8	0.26	166.3	2.00	0.2500	0.010	2.8
	L0012N0376	14.4	3763	5600.0	16.84	2.1$\times {10}^3$	64.01	16.0000	0.025	2.8
	L0012N0376	14.4	3763	2800.0	8.42	10.6$\times {10}^3$	32.01	8.0000	0.025	2.8
	L0012N0376	14.4	3763	1400.0	4.21	5320.0	16.00	4.0000	0.025	2.8
	L0012N0376	14.4	3763	700.0	2.11	2660.0	8.00	2.0000	0.025	2.8
	L0012N0376	14.4	3763	350.0	1.05	1330.0	4.00	1.0000	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	175.0	0.53	665.0	2.00	0.5000	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	87.5	0.26	332.5	1.00	0.2500	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	43.8	0.13	166.3	0.50	0.1250	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	21.9	0.07	83.13	0.25	0.0625	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	10.9	0.03	41.6	0.13	0.0313	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	5.5	0.02	20.8	0.06	0.0156	0.025, 0.0025	2.8
	L0012N0376	14.4	3763	350.0	1.05	350.0	1.05	1.0000	0.025	0.0
	L0012N0376	14.4	3763	350.0	1.05	2609.2	8.09	1.0000	0.025	10.0
	L0012N0188	115.0	1883	700.0	1.05	2660.0	4.00	1.0000	0.025	2.8
	L0012N0188	115.0	1883	350.0	0.53	1330.0	2.00	0.5000	0.025	2.8
	L0012N0188	115.0	1883	175.0	0.26	665.0	1.00	0.2500	0.025	2.8
	L0012N0188	115.0	1883	87.5	0.13	332.5	0.50	0.1250	0.025	2.8
	L0012N0188	115.0	1883	43.8	0.07	166.3	0.20	0.0625	0.025	2.8
	L0012N0188	115.0	1883	21.9	0.03	83.13	0.13	0.0313	0.025	2.8

Equations

$$\rho (r)=\frac{{\rho }_s}{r/r_s(1+r/r_s)},$$

$$l\equiv {\mathrm{L}}_{\mathrm{b}}/N_{\mathrm{p}}^{1/3},$$

$${ϵ}_{\mathrm{sp}}\equiv 2.8\times ϵ.$$

$$M_{\lim }\approx 10.1\overline{\rho }l{k}_{\mathrm{peak}}^2,$$

$$\begin{split}r_{\rm lim}&\approx 1.3\,l\,\biggr{(}\frac{\Omega_{\rm DM}}{\Delta}\biggl{)}^{1/3}\biggr{(}\frac{k_{\rm Ny}}{k_{\rm peak}}\biggl{)}^{2/3},\\ &\approx 0.07\,l\,\biggr{(}\frac{k_{\rm Ny}}{k_{\rm peak}}\biggl{)}^{2/3},\end{split}$$

$$\begin{split}r_{200}=&\biggl{(}\frac{3\,\Omega_{\rm DM}}{8\,\pi}\biggr{)}^{1/3}\,\biggl{(}\frac{N_{200}}{100}\biggr{)}^{1/3}\,l\\ \approx&\,\,0.332\,l\,\biggl{(}\frac{N_{200}}{100}\biggr{)}^{1/3}.\end{split}$$

$$\begin{split}\epsilon_{v}&\gg l\,\biggl{(}\frac{3\,\Omega_{\rm DM}}{800\,\pi}\frac{1}{N_{200}^{2}}\biggr{)}^{1/3}\\ &\approx 3.32\times 10^{-3}\,l\,\biggr{(}\frac{N_{200}}{100}\biggl{)}^{-2/3}.\end{split}$$

$$|\delta v_{\perp }|\approx \frac{2G{m}_{\mathrm{DM}}b}{(b^2+{ϵ}^2)v},$$

$$\begin{split}\Delta v_{\perp}^{2}&=8\,N\biggr{(}\frac{Gm_{\rm DM}}{Rv}\biggl{)}^{2}\int_{b_{\rm min}}^{b_{\rm max}}b^{3}(b^{2}+\epsilon^{2})^{-2}db\\ &=\frac{4v^{2}}{N}\biggr{[}\ln(\epsilon^{2}+b^{2})+\frac{\epsilon^{2}}{\epsilon^{2}+b^{2}}\biggl{]}_{b_{\rm min}}^{b_{\rm max}},\end{split}$$

$$\frac{t_{\mathrm{relax}}}{t_{\mathrm{orb}}}=\frac{N}{8\ln (b_{\max }/b_{\min })}\simeq \frac{N}{8\ln (N/2)},$$

$$\frac{t_{\mathrm{relax}}}{t_{\mathrm{orb}}}$$

$$\displaystyle\approx\frac{N}{4}\biggr{[}\ln\biggr{(}\frac{R^{2}}{\epsilon^{2}}\biggl{)}-\ln\biggr{(}\frac{3}{2}\biggl{)}-\frac{2}{3}\biggl{]}^{-1}$$

$$\approx \frac{N}{8\ln (R/ϵ)},$$

$$\begin{split}\kappa_{\rm P03}(r,z)&\equiv\frac{t_{\rm relax}(r)}{t_{\rm 200}(z)}\\ &=\frac{N}{8\ln N}\frac{r/V_{c}(r)}{r_{200}(z)/V_{200}(z)}\\ &=\frac{\sqrt{200}}{8}\frac{N}{\ln N}\biggr{(}\frac{\overline{\rho}(r)}{\rho_{\rm crit}(z)}\biggl{)}^{-1/2}\end{split}$$

$$\kappa_{\epsilon}(r,z)=\frac{\sqrt{200}\,N}{4}\,\biggr{[}\ln\biggr{(}\frac{R^{2}}{\epsilon^{2}}+1\biggl{)}+\\ \frac{\epsilon^{2}-2R^{2}}{3(\epsilon^{2}+R^{2})}-\ln\biggr{(}\frac{2}{3}\biggl{)}\biggl{]}^{-1}\biggr{(}\frac{\overline{\rho}(r)}{\rho_{\rm crit}(z)}\biggl{)}^{-1/2}.$$

$$\frac{r_{\mathrm{conv}}}{r_{200}}=\frac{{\kappa }_{P03}^{2/3}C}{N_{200}^{1/3}},$$

$$r_{\mathrm{conv}}(z)$$

$$\displaystyle=C\,\kappa_{\rm P03}^{2/3}\,\biggr{(}\frac{3\,\Omega_{\rm DM}}{800\,\pi}\biggl{)}^{1/3}\,l(z),$$

$$r_{\mathrm{conv}}$$

$$\approx 5.5\times 10^{-2}l(z),$$

$$\frac{\Delta V_c}{V_c}(\psi )=-\log (1-10^{a{\psi }^2+b\psi +c}),$$

Full text

11footnotetext: E-mail: aaron.ludlow@icrar.org

1 Introduction

Cosmological simulations have become an essential component of astronomical science. Simulations of collisionless cold dark matter (CDM), in particular, have matured to a point where both the statistical properties of large-scale structure, as well as the highly non-linear structure of dark matter haloes are largely agreed upon, even between groups employing widely varying simulation or analysis methods. Among these are: the topology of large-scale structure (e.g. Gott et al., 1987; James et al., 2007; Blake et al., 2014); the matter power spectrum (e.g. Smith et al., 2003); the clustering (e.g. Kaiser, 1984; White et al., 1987; Poole et al., 2015; Tinker et al., 2010), mass function (e.g. Jenkins et al., 2001; Reed et al., 2003; Tinker et al., 2008; Despali et al., 2016) and shapes (e.g. Allgood et al., 2006; Despali et al., 2014; Vera-Ciro et al., 2014; Vega-Ferrero et al., 2017) of dark matter haloes; their spherically averaged mass profiles (e.g. Navarro et al., 1996, 1997; Bullock et al., 2001; Ludlow et al., 2013; Dutton & Macciò, 2014); mass assembly histories (e.g. van den Bosch, 2002; Zhao et al., 2009; Correa et al., 2015a, b) and the mass function and radial distribution of their substructure populations (e.g. Ghigna et al., 1998; Stoehr et al., 2003; Gao et al., 2004; Springel et al., 2008).

The radial mass profile of dark matter haloes is a particularly important and robust prediction of N-body simulations. For relaxed haloes it can be approximated by the NFW profile (Navarro et al., 1996, 1997), though slight deviations from this form have been reported extensively in the literature (e.g. Navarro et al., 2004; Ludlow et al., 2013; Dutton & Macciò, 2014; Ludlow & Angulo, 2017). The NFW profile has a central cusp where densities diverge as $\rho\propto r^{-1}$ and a steep outer profile where $\rho(r)$ tapers off as $r^{-3}$. In parametric form, the spherically averaged density profiles can be well approximated by

[TABLE]

where $\rho_{s}$ and $r_{s}$ are characteristic values of density and radius.

Agreement on these issues required pain-staking tests of numerical convergence that demanded repeatability of simulation results, regardless of the numerical methods employed or the numerical parameters adopted. A number of these studies led to the development of useful “convergence criteria” that can be used to disentangle aspects of simulations that are reliably modelled from those that may be affected by numerical artifact. These studies differ in the details, but uniformly agree that systematic convergence tests are necessary to validate the robustness of a particular numerical result. Numerical requirements for convergence in halo mass functions, for example, may differ substantially from those required for convergence in shapes, dynamics or mass profiles.

For collisionless CDM, once a cosmological model has been specified, only numerical parameters remain. The starting redshift and finite box size of simulations, for example, affect halo mass functions and clustering, but leave the internal properties of dark matter haloes largely unchanged (Knebe et al., 2009; Power & Knebe, 2006). Other parameters impose strict limits on spatial resolution, or otherwise alter the inner structure of haloes in non-trivial ways. Of particular importance are the gravitational force softening, $\epsilon$ (which prevents divergent pairwise forces and suppresses large-angle deflections), the integration timestep for the equations of motion, $\Delta t$, and the particle mass resolution, $m_{p}$.

Power et al. (2003, hereafter P03) provided a comprehensive survey of how these numerical parameters affect the internal structure of a simulated CDM halo. They conclude that convergence in mass profiles can be achieved for suitable choices of timestep, softening, particle number and force accuracy. For choices of softening that suppress discreteness effects, and for timesteps substantially shorter than the local dynamical time, circular velocity profiles converge to within \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10 per cent roughly at the radius enclosing a sufficient number of particles to ensure that the local 2-body relaxation time exceeds a Hubble time. Their tests led to the development of what is now a standard choice for the “optimal” softening for cosmological simulations, and to an empirical prescription for calculating the “convergence radius” of dark matter haloes. Their results – which we put to the test in subsequent sections – have been validated and extended by a number of follow-up studies (e.g. Diemand et al., 2004; Springel et al., 2008; Navarro et al., 2010; Gao et al., 2012).

Strictly speaking, the criteria laid out by P03 mainly apply to convergence in the circular velocity profiles, $V_{c}(r)$, of individual haloes, and may not apply to convergence in other quantities of interest, such as their shapes (Vera-Ciro et al., 2014), mass functions (e.g. Reed et al., 2003), to various aspects of their substructure distributions (e.g. Ghigna et al., 2000; Reed et al., 2005; Springel et al., 2008), or to population-averaged profiles of, for example, density or circular velocity. We focus on the latter in this paper.

P03 defined convergence empirically: the same simulation was repeated multiple times using different numerical parameters and the results were used to quantify the radial range over which $V_{c}(r)$ remained insensitive to those choices. van den Bosch & Ogiya (2018) follow a different approach. Using a series of idealized numerical experiments, they argue that inappropriate choices for gravitational softening are detrimental to the evolution of substructure haloes and that, as a result, many state-of-the-art cosmological simulations are still subject to the classic “over-merging” problem (Moore et al., 1996). They additionally argue that discreteness-driven instabilities in subhaloes with \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{3} particles forbids a proper assessment of their evolution in strong tidal fields, limiting our ability to interpret convergence in their mass functions or internal structure.

In this paper, we interpret convergence in the median $V_{c}(r)$ profiles of halos in terms of 2-body relaxation, which varies from system-to-system depending on the number of particles used to sample their mass profile. To do so, we carry out the same simulation several times but with different numbers of particles, allowing us to compare the same sub-population (i.e. those occupying a particular mass bin) at varying resolution. It is important to note, however, that 2-body scattering may not be the sole driver of relaxation in collisionless systems. Collective relaxation, for example, is a distinct process by which individual particle trajectories are altered repeatedly in response to large-scale and potentially long-lived fluctuations to the global potential (see e.g. Weinberg, 1993). This is neglected in simple analytic estimates of collisional relaxation rates, such as those presented in Section 3.2 (see also Chandrasekhar, 1943; Hénon, 1961; Cohn & Kulsrud, 1978), but may dominate when fluctuations are comparable to the size of the system under consideration. The fluctuations may be physical, or form as a result of discreteness-driven noise in collisionless $N$-body simulations. Indeed, it has been argued that this process may give rise to artificial fragmentation that plagues traditional $N$-body simulations of warm dark matter models, but appears to have measurable effects in CDM simulations as well (see, e.g., Power et al., 2016, and references therein).

Going beyond pure dark matter (DM), hydrodynamic simulations of galaxy formation are reaching new levels of maturity. The increase of computational resources and improved algorithms enable fully cosmological simulations of galaxy formation to be carried out; simulations that often resolve tens of thousands of individual galaxies in volumes approaching those required for cosmological studies. Notable among these are eagle (Schaye et al, 2015; Crain et al., 2015), Illustris (Vogelsberger et al., 2014a; Pillepich et al, 2018), Horizon-AGN (Dubois et al, 2014), and the Magneticum Pathfinder simulations (Dolag et al., 2016). Although sub-grid models in simulations such as these must be calibrated to reproduce a desired set of observables (for eagle, the $z=0$ galaxy stellar mass function and size-mass relation), many of their predictions have been ratified by observations making them a useful tool for interpreting observational data and for illuminating the complex physical processes that give rise to galaxy scaling relations.

As discussed by Schaye et al (2015), the need to calibrate subgrid models in simulations affects our ability to interpret numerical convergence, particularly when mass and spatial resolution are improved. Clearly convergence is a requirement for predictive power, but for a multi-scale process such as galaxy formation it should arguably be attained only after recalibration of the subgrid physics, thus allowing models to benefit from increased resolution by incorporating new, scale-dependent physical processes.

The convergence of hydrodynamic simulations is, in any event, poorly understood. It remains unclear, for example, how robust predictions of galaxy formation models are to small changes in the numerical parameters. To cite an example, Ludlow et al. (2019) recently pointed out that simulated galaxy sizes, quantified by their projected half-mass radii, grow over time as a result of spurious energy transfer between stellar and DM particles of unequal mass. The effect can be suppressed by adopting a stellar to DM particle mass ratio that is close to unity, even when other numerical and subgrid parameters are unchanged. We seek to further address these issues using a suite of simulations drawn from the eagle project. In this first paper we study the sensitivity of our simulations to numerical parameters when only dark matter is present, seeking to illuminate and clarify shortcomings of prior work. In a follow-up paper, we will address convergence in fully-hydrodynamical simulations using a well-calibrated galaxy formation model from the eagle project.

Our study is part of the eagle Project. All of our runs were carried out using the same simulation code and adopt the same “fiducial” numerical parameters as Schaye et al (2015), which we vary systematically from run-to-run. We concentrate our discussion primarily on gravitational softening and the impact of 2-body collisions on the spatial resolution of N-body simulations, but consider mass resolution and timestepping when necessary. Softening has been studied in great detail in the past two decades, but this has not led to a clear consensus on what constitutes an “optimal” softening length for a given simulation.

The remainder of the paper is organized as follows. In Section 2.1 we describe our simulation suite and the variation of numerical parameters, as well as halo finding techniques (Section 2.2). We provide simple analytic estimates of plausible bounds on gravitational softening in Section 3; we introduce the 2-body “convergence radius” in Section 3.2, highlighting its explicit dependence on softening. We then present our main results: the convergence of median circular velocity profiles is discussed in Section 4, followed by that of mass functions (Section 5). We summarize and conclude in Section 6.

2 Simulations

2.1 Simulation set-up

All runs were carried out in the same ${\rm L_{b}}=12.5\,{\rm Mpc}$ cubic periodic volume which was simulated repeatedly using different numbers of particles, $N_{\rm p}$, gravitational softening lengths, $\epsilon$, and timestep size, $\Delta t$. Each run adopted the set of “Planck” cosmological parameters used for eagle (Schaye et al, 2015; Planck Collaboration et al., 2014): $\Omega_{\rm M}=\Omega_{\rm DM}+\Omega_{\rm bar}=1-\Omega_{\Lambda}=0.307$; $\Omega_{\rm bar}=0.04825$; $h=0.6777$; $\sigma_{8}=0.8288$; $n_{s}=0.9611$. Here $\Omega_{i}$ is the energy density of component $i$ expressed in units of the critical density, $\rho_{\rm crit}\equiv 3H_{0}^{2}/(8\pi G)$; $h\equiv H_{0}/[100\,{\rm km/}s/{\rm Mpc}]$ is Hubble’s constant; $\sigma_{8}$ is the $z=0$ linear rms density fluctuation in 8 $h^{-1}{\rm Mpc}$; and $n_{s}$ is the primordial power spectral index. Initial conditions for each simulation were generated using second-order Lagrangian perturbation theory at $z=127$ (Jenkins, 2013), which is sufficiently high to ensure that all resolved modes are initially well within the linear regime. We sample the linear density field with $N_{\rm p}=188^{3}$, $376^{3}$ and $752^{3}$ equal-mass particles; the corresponding particle masses are $m_{\rm DM}=1.15\times 10^{7}$, $1.44\times 10^{6}$ and $1.80\times 10^{5}\,{\rm M}_{\odot}$. All simulations sample resolved modes using the same initial phases. For consistency we label each run in a way that encodes the box size and particle number, using the same nomenclature as Schaye et al (2015). For example, L0012N0376 corresponds to a run with ${\rm L_{b}}=12.5\,{\rm Mpc}$ and $N_{\rm p}=376^{3}$ particles. The DM density field is evolved using the same version of P-gadget (Springel, 2005) employed for the eagle project.

It is common in the literature to refer to $\epsilon$ as the “spatial resolution” of a simulation, not surprisingly given its dimensions. For cosmological simulations of uniform mass resolution it is customary to adopt a gravitational softening length that is a fixed fraction of the (comoving) mean inter-particle separation,

[TABLE]

thus fixing the ratio $\epsilon/m_{\rm DM}^{1/3}$. In eagle the softening parameter, initially fixed in comoving coordinates, reaches a maximum physical value at redshift $z_{\rm phys}=2.8$, after which it remains constant in physical coordinates. For ${\rm L_{b}=100}\,{\rm Mpc}$, $N_{\rm p}=1504^{3}$ and $\epsilon(z=0)=700\,{\rm pc}$, this implies, $\epsilon_{\rm phys}/l\approx 0.011$ for $z\leq 2.8$, while $\epsilon_{\rm cm}/l\approx 0.04$ in co-moving coordinates for $z>2.8$. We will hereafter refer to $\epsilon_{\rm fid}(z=0)/l\approx 0.011$ as the “fiducial” softening length regardless of mass resolution, and will vary $\epsilon$ away from this value by successive factors of two. For ${\rm N_{\rm p}=752^{3}}$ the “fiducial” softening length is $\epsilon_{\rm fid}=175\,{\rm pc}$; $\epsilon_{\rm fid}=350\,{\rm pc}$ for $N_{\rm p}=376^{3}$, and $700\,{\rm pc}$ for $N_{\rm p}=188^{3}$. We further note that $\epsilon$ refers to the Plummer-equivalent softening length, which is related to the “spline” softening length used by P-gadget through

[TABLE]

To test the sensitivity of our results to changing $z_{\rm phys}$, we have also carried out runs with $z_{\rm phys}=0$ (fixed co-moving $\epsilon$ at all $z$) and $\infty$ (fixed physical $\epsilon$). For convenience, we will sometimes reference softening parameters relative to the fiducial value, hence defining the relative softening length $f_{\epsilon}\equiv\epsilon_{i}/\epsilon_{\rm fid}$. Table 1 lists all of the relevant numerical parameters for our simulations.

2.2 Halo identification

We identify haloes in all of our simulations using the subfind (Springel et al., 2001) algorithm. subfind first links dark matter particles into friends-of-friends (FoF) groups before separating them into a number of self-bound “subhaloes”. Each FoF group contains a central or “main” subhalo that contains most of its mass, and a number of lower-mass substructures. For each FoF halo and its entire hierarchy of subhaloes subfind records a number of attributes, the most basic of which include its mass, ${\rm M_{FoF}}$ (for FoF haloes), position (defined as the location of the particle with the minimum potential energy) and the magnitude and location of its peak circular speed, $V_{\rm max}$ and $r_{\rm max}$. For FoF haloes (defined as “main” haloes in what follows), subfind also records other common mass definitions based on a variety of spherical overdensity (SO) boundaries: ${\rm M_{200}}$ is the mass contained within a spherical region whose mean density is 200$\times\rho_{\rm crit}(z)$; ${\rm M_{\Delta}}$ encloses a mean density of $\Delta\times\rho_{\rm crit}(z)$, where $\Delta$ is the redshift-dependent virial overdensity of Bryan & Norman (1998) ($\Delta=103.7$ for our adopted cosmology). Using as a starting-point all FoF particles in a given halo, subfind also computes its gravitationally-bound mass, ${\rm M_{bound}}$, as well as the mass bound to each of its subhaloes. Each of these are common and sensible ways to define halo masses, which we compare in Section 5.2. Note that the virial mass of a halo implicitly defines its virial radius, $r_{\Delta}$, and corresponding circular velocity, ${\rm V_{\Delta}}$: for an overdensity $\Delta$, for example, $r_{\Delta}=3M_{\Delta}/4\pi\Delta\rho_{\rm crit}$, and $V_{\Delta}=\sqrt{G\,M_{\Delta}/r_{\Delta}}$.

In addition to halo properties, subfind also records lists of particles belonging to each halo (which include all bound as well as unbound particles within $r_{200}$), which we use to calculate their spherically-averaged enclosed mass profiles, $M(r)$. Note that all particles contribute to $M(r)$, and not only those that subfind deems bound to the halo. This choice precludes any subtle dependence of our results on subfind’s unbinding algorithm. We construct these profiles for all main haloes in 50 equally-spaced logarithmic bins spanning $-5\leq\log r/[{\rm Mpc}]\leq 0$ which we then use to build median circular velocity profiles, ${\rm V_{c}}(r)=\sqrt{G\,M(r)/r}$, as a function of halo mass, and various other structural scaling relations. Note that all particles are used to calculate $M(r)$, and not just those deemed bound to the main halo by subfind. The results presented in the following sections are obtained using these spherically-averaged profiles, and the halo masses returned by subfind.

Assessing convergence in properties of substructure is challenging (see van den Bosch & Ogiya, 2018; van den Bosch et al., 2018, for a recent in depth analysis). For that reason, we will focus our analysis on isolated haloes expected to host central galaxies in hydrodynamic runs.

3 Analytic expectations

3.1 Preliminaries: limits on gravitational softening

Softening gravitational forces in N-body simulations has distinct advantages. In particular, it suppresses large-angle deflections due to (artificial) 2-body scattering, thereby permitting the use of low-order orbit integration schemes. This substantially decreases wall-clock times required for a given N-body problem. There are, however, drawbacks: when $\epsilon$ is small shot noise in the particle load can result in large near-neighbor forces, or when large, to systematic suppression of inter-particle forces. Both effects can jeopardize the results of N-body simulations and the optimal choice of gravitational softening should provide a compromise between the two.

The finite particle mass and limited force resolution inherent to particle-based simulations can give rise to adverse discreteness effects if $\epsilon$ is not properly chosen, and the debate over what constitutes a wise choice persists. Some studies, designed specifically to annotate discreteness-driven noise in simulations, suggest a safe lower limit to softening of \epsilon/l\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}1-2 (e.g. Melott et al., 1997; Power et al., 2016). This is supported by others who argue that various two-point statistics of the DM density field are not converged on scales \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}l (Splinter et al., 1998). It is important to note, however, that discreteness noise does not propagate from the small scales where it is introduced to larger ones, being typically confined to scales of order $\sim\epsilon$ to $\sim 2\,l$ (Romeo et al., 2008).

Whether simulations are trustworthy below scales $\sim l$ remains a matter of debate. Claims above to the contrary are often based on particular cases: cold, plane-symmetric collapse and simulations with truncated linear power spectra give rise to spurious clustering as a result of discreteness-driven relaxation. The effect is prominent in simulations of hot and warm DM models (e.g. Avila-Reese et al., 2001; Bode et al., 2001; Knebe et al., 2003), where counterfeit haloes form as a result of anisotropic force errors occurring during the initial phases of filamentary collapse (Hahn et al., 2013). These haloes follow a spatial pattern that reflects the Lagrangian inter-particle spacing, but remain prominent in warm DM simulations that adopt isotropic softening lengths of order $l$ (Power et al., 2016). Better success at stifling these haloes has been achieved using adaptive, anisotropic softening lengths (Hobbs et al., 2016) or alternatives to traditional N-body methods (Angulo et al., 2013).

Power et al. (2016) argue that spurious haloes are also present in CDM models, but their prevalence has not been quantified due to difficulties disentangling them from “genuine” ones. Using simulations of hot DM models, Wang & White (2007) showed that spurious structures dominate the halo mass function on scales below a limiting mass given by

[TABLE]

where $\overline{\rho}=\Omega_{\rm DM}\,\rho_{\rm crit}$ is the mean matter density and $k_{\rm peak}$ is the wavenumber at which the dimensionless matter power spectrum, $k^{3}\,P(k)$, reaches a maximum. Assuming $M_{\rm lim}=(4/3)\,\pi\,r_{\rm lim}^{3}\,\Omega_{\rm DM}\Delta\rho_{\rm crit}$, we can write eq. 4 in terms of the maximum comoving size of spurious haloes:

[TABLE]

where $k_{\rm Ny}\equiv\pi/l$ is the Nyquist frequency, and in the second expression we have used $\Omega_{\rm DM}=0.307$ and $\Delta=200$. For CDM models, finite resolution suggests that $k_{\rm peak}\approx k_{\rm Ny}$, implying $r_{\rm lim}\approx 0.07\,l$ (assuming eq. 4 remains valid in this regime).

It is useful to compare $r_{\rm lim}$ to the comoving virial radius of the lowest-mass haloes resolved in cosmological CDM simulations. Noting that $M_{200}=N_{200}\,m_{\rm DM}$ and $m_{\rm DM}=\rho_{\rm crit}\,\Omega_{\rm DM}\,l^{3}$, this can be written

[TABLE]

Note, for CDM, $r_{200}=r_{\rm lim}\approx 0.07\,l$ when $N_{200}\approx 1$, suggesting spurious haloes are negligible for models in which $k^{3}\,P(k)$ increases to arbitrarily small scales. For models with suppressed small-scale power, on the other hand, spurious haloes may dominate at low masses. For example, r_{\rm lim}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}l provided the Power spectrum “turn over” is resolved with \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}54 resolution elements.

To avoid biasing gravitational collapse at the resolution limit of CDM simulations, the comoving softening length should therefore remain smaller than the comoving virial radius of the lowest-mass haloes that can be resolved. Indeed, Lukić et al. (2007) and Power et al. (2016) showed that CDM halo mass functions are substantially suppressed on scales below which the softening length exceeds the halo virial radius. For a conservative resolution limit of $N_{200}=100$, eq. 6 suggests that $\epsilon$ should remain smaller than about one third of the mean inter-particle spacing; for $N_{200}=20$, \epsilon/l\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.2 is required.

Most recent large-scale simulations adopt softening lengths substantially smaller than these limits, but large enough to ensure that the lowest-mass haloes resolved by the simulation remain approximately collisionless***No haloes are truly collisionless in particle-based simulations. We use the phrase here to refer to haloes in which collisional dynamics driven by 2-body interactions remain small compared to those dictated by the global potential. at all times. One such criterion demands that the specific binding energy of low-mass haloes, $V_{200}^{2}\simeq(10\,G\,H\,N_{200}\,m_{\rm DM})^{2/3}$, remains larger than the binding energy of two DM particles separated by $\epsilon$: $v_{\epsilon}^{2}=G\,m_{\rm DM}/\epsilon$. The condition $v_{\epsilon}^{2}\ll V_{200}^{2}$ imposes a lower limit on $\epsilon$ of

[TABLE]

Eq. 7 can also be expressed $\epsilon_{v}>r_{200}/N_{200}$, where we have used the fact that $M_{200}=\Omega_{\rm DM}\,\rho_{\rm crit}\,N_{200}\,l^{3}$. This is comparable to the softening length required to suppress large-angle deflections during close encounters, given by $\epsilon_{90}\sim b_{90}=2Gm_{\rm DM}/v^{2}$, where $b_{90}$ is the impact parameter giving rise to $\sim 90^{\circ}$ deflections (Binney & Tremaine, 1987) and $v^{2}\approx GM/r$ is the typical speed of particles at distance $r$. This condition therefore requires \epsilon_{90}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}2\,r_{200}/N_{200}, a factor of 2 larger than $\epsilon_{v}$. Suppressing large accelerations during to close encounters results in stricter limits on softening (see P03). For example, requiring that the maximum stochastic acceleration due to close encounters, $a_{\epsilon}=Gm_{\rm DM}/\epsilon^{2}$, remains smaller than the minimum mean-field acceleration across the entire halo, $a_{\rm min}=G\,M_{200}/r_{200}^{2}$, imposes a lower limit of \epsilon_{\rm acc}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}r_{200}/\sqrt{N_{200}}.

For $N_{vir}\sim 20$ eq. 7 implies a lower limit of \epsilon/l\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.01, comparable to values adopted for essentially all recent large-scale simulations projects. The Bolshoi simulation (Klypin et al., 2011), for example, had a force resolution of $\approx 0.016\,l$, while the Multi-dark simulations (Klypin et al., 2016) adopted $\epsilon/l=0.014-0.026$ (Plummer equivalent); the Millennium (Springel et al., 2005), Millennium-II (Boylan-Kolchin et al., 2009) and Millennium-XXL (Angulo et al., 2012) each used $\epsilon/l\approx 0.022$; $\epsilon/l=0.016$ for DUES FUR (Alimi et al., 2012); $\epsilon/l=0.020$ for Copernicus Complexio (Hellwing et al., 2016). Cosmological hydrodynamical simulations use comparable values of softening: as mentioned above, eagle adopted a physical softening length of $\epsilon/l=0.011$ ($\epsilon/b=0.04$ in co-moving coordinates at $z>2.8$), while Illustris-TNG used maximum physical value of $\epsilon/l\approx 0.012$ (Springel et al., 2018). All values above are quoted as physical softening lengths at $z=0$, unless stated otherwise. We will see in Section 3.2 that the median circular velocity profiles of haloes in collisionless cosmological simulations are affected by 2-body scattering at radii that generally exceed these softening lengths. For example, $V_{c}(r)$ converges to better than 10 per cent a radii r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.055\times l, and to better than 3 per cent for r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.10\times l (section 4.3). In all of these runs, softening is unlikely to compromise the mass profiles of haloes on scales not already influenced by 2-body relaxation.

As mentioned above, the “optimal” softening, $\epsilon_{\rm opt}$, for a given simulation is the one that balances large force errors due to shot noise with biases resulting from large departures from the Newtonian force law. Merritt (1996) suggested that $\epsilon_{\rm opt}$ should be chosen to minimize the average square error in force evaluations relative to what is expected from an equivalent smooth matter distribution. A drawback of this approach is that $\epsilon_{\rm opt}$ depends not only on the mass distribution under consideration – which is not generally known a priori – but also on the number of particles in the system, $N$. Dehnen (2001), for example, found that the optimal softening for a Hernquist (1990) halo is roughly $\epsilon_{\rm opt}/a\simeq 0.017\,(N/10^{5})^{-0.23}$, where $a$ is its scale radius. van den Bosch & Ogiya (2018) find that $\epsilon_{\rm opt}/r_{200}=0.005\times(N_{200}/10^{5})^{-1/3}$, a result obtained by studying the stability of the central cusp of an isolated NFW halo ($N_{200}=10^{5}$ and $c=r_{200}/r_{s}=10$) to long-term secular evolution ($t\approx 60\,{\rm Gyr}$); the scaling with $N_{200}$ is motivated by the work of van Kampen (2000). We will see in Section 5.3 that this is sufficiently small to avoid compromising the spatial resolution in halo centres in cosmological simulations using equal-mass particles. In our study, we consider a broad range of softening lengths, spanning $\epsilon/l\approx 0.17$ to $\approx 1.6\times 10^{-4}$.

Figure 1 compares the softening lengths used in a number of recent state-of-the-art simulations, and shows for comparison the characteristic values $\epsilon_{v}$, $\epsilon_{90}$ and $\epsilon_{acc}$ mentioned above. Note that, in all cases, $\epsilon$ is expressed in units of the virial radius, $r_{200}$, and plotted as a function of $N_{200}$.

In summary, a number of previous studies suggest that there can be considerable errors affecting the results of N-body simulations on scales \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}l (e.g. Peebles et al., 1989; Melott et al., 1997; Splinter et al., 1998; Romeo et al., 2008). These may be driven by discreteness-noise (e.g. Power et al., 2016) or errors that affect phases of the Fourier component on such scales (see, e.g., Melott, 2007). Their severity, however, depends on the linear power spectrum adopted, or the type of simulation or statistic studied. Regardless, this underscores the need to establish independent convergence criteria for all results obtained from simulations, and to abstain from the common view that the softening length, or particle mass, represent meaningful “resolution limits”.

3.2 2-Body relaxation and the convergence radius

3.2.1 The relaxation timescale

It is important to emphasize that softening pair-wise forces does not necessarily increase 2-body relaxation times, which generally impose stricter constraints on the spatial resolution of N-body simulations. To see why, consider the cumulative effect of 2-body interactions incurred by a test particle as it crosses an $N$-particle system with surface mass density $\Sigma\approx N/\pi R^{2}$. Following Binney & Tremaine (1987) (see also Huang et al., 1993; Farouki & Salpeter, 1982), we assume that any one encounter induces a small velocity perturbation $|\delta v_{\perp}|\ll v$ perpendicular to the particle’s direction of motion, but leaves its trajectory unchanged. The perturbation due to a single encounter can be expressed

[TABLE]

where $b$ is the impact parameter and $\epsilon$ the (Plummer) softening length. In a single crossing, the test particle will experience $\delta n\approx 2\pi\Sigma\,b\,db$ such collisions with impact parameters spanning $b$ to $b+db$. Integrating the square velocity change over all such encounters yields

[TABLE]

where $b_{\rm min}$ and $b_{\rm max}$ are, respectively, the minimum and maximum impact parameters, and we have assumed a typical velocity $v^{2}=G\,m_{\rm DM}\,N/R$.

The relaxation time can be defined in terms of the number of orbits a test particle must execute in order to lose memory of its initial trajectory, i.e. $n_{\rm orb}\sim v^{2}/\Delta v_{\perp}^{2}\sim t_{\rm relax}/t_{\rm orb}$. In the limit $b\gg\epsilon$, eq. 9 implies

[TABLE]

where $t_{\rm orb}=2\,\pi\,r/v$ is the local orbital time, and we have assumed $b_{\rm max}=R$ and $b_{\rm min}=b_{90}\equiv 2G\,m_{\rm DM}/v^{2}$. Eq. 10 is comparable to the classic relaxation time calculated by Binney & Tremaine (1987) for Keplerian forces, and depends only on the enclosed number of particles. Note that for a narrow range of impact parameter, $b_{\rm max}/b_{\rm min}=(1+\Delta)$, eq. 10 implies $t_{\rm relax}/t_{\rm orb}\propto\ln(1+\Delta)^{-1}$: fixed intervals of $b$ therefore contribute equally to relaxation, which is sensitive to both close and distant encounters. As a result, softening forces on small scales does little to prolong relaxation, which is primarily driven by large numbers of distant perturbers.

What is the appropriate choice of $b_{\rm min}$ for haloes in cosmological simulations? If $\epsilon$ is chosen in order to suppress large-angle deflections, then the assumption $b_{\rm min}=b_{90}$ may need revision. Indeed, for a Plummer potential there is a well-defined impact parameter†††This can be seen by maximizing the radial gradient of a Plummer force law and solving for $\epsilon$. The minimum impact parameter obtained this way is slightly smaller than that which maximizes the velocity perturbation in eq. 8, which occurs at $b=\epsilon$., $b_{\rm min}=\epsilon/\sqrt{2}$, for which the perpendicular pair-wise force between particles is maximized, tending to zero for both larger and smaller $b$. Inserting this into eq. 9, and assuming $b_{\rm max}=R$, we obtain

[TABLE]

where the last two steps are valid provided $R\gg\epsilon$. Note that eq. 13 depends logarithmically on $\epsilon$, and is equivalent to eq. 10 if $b_{\rm max}=R$ and $b_{\min}=\epsilon$. This confirms our expectation that softened forces give rise to modest changes in $t_{\rm relax}$, even for very small values of the softening length (see also, e.g., Huang et al., 1993; Theis, 1998; Dehnen, 2001; Diemand et al., 2004). We will test this explicitly in Section 4.3.

3.2.2 The convergence radius

Collisional relaxation imposes a lower limit on the spatial resolution of N-body simulations that is typically larger than the limits on gravitational softening outlined in the previous section. Based on an extensive suite of simulations of individual Milky Way-mass haloes, P03 concluded that convergence is obtained at radii that enclose a sufficient number of particles so that the local two-body relaxation timescale is comparable to or longer than a Hubble time, $t_{H}(z)\sim t_{200}(z)$, where $t_{200}(z)=2\,\pi\,r_{200}(z)/V_{200}(z)$ is the circular orbital time at $r_{200}$. The “convergence radius”, $r_{\rm conv}$, implied by eq. 10 can thus be approximated by the solution to

[TABLE]

where $\overline{\rho}(r)$ is the mean internal density enclosing $N\equiv N(r)$ particles. We can rewrite eq. 14 in order to incorporate its explicit dependence on softening implied by eq. 11, resulting in

[TABLE]

Note that we have used subscripts on eqs. 14 and 15 in order to distinguish the P03 convergence radius from our softening-dependent estimate above, a convention we retain throughout the paper. Either equation, however, can be used to estimate $r_{\rm conv}$ once an empirical relationship between $\kappa\equiv t_{\rm relax}/t_{\rm H}$ and “convergence” has been determined from simulations. P03 found that circular velocity profiles converge to better than $10$ per cent provided \kappa_{\rm P03}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.6. Navarro et al. (2010, hereafter N10) showed that better convergence can be obtained for larger values of $\kappa_{\rm P03}$; at $\kappa_{\rm P03}=7$, for example, ${\rm V_{c}}(r)$ profiles converge to $\approx 2.5$ per cent. We follow P03 and N10 and quantify convergence in circular velocity profiles using $\Delta V_{c}\equiv(V_{c}^{\rm high}-V_{c}^{\rm low})/V_{c}^{\rm low}$, where $V_{c}^{\rm low}$ and $V_{c}^{\rm high}$ are the median profiles for haloes of fixed mass in our low- and highest-resolution simulations, respectively.

3.2.3 A simpler convergence criterion

As a useful approximation, we can rewrite the convergence radius implied by eq. 14 as

[TABLE]

where $C\equiv 4\,(\ln N_{c}/\sqrt{N_{c}})^{2/3}$ and $N_{c}\equiv N(<r_{\rm conv})$; these quantities depend weakly on concentration and $N_{200}$, but also on $\kappa_{\rm P03}$. To illustrate, we set $\kappa_{\rm P03}=1$ and solve eq. 14 assuming an NFW mass profile to determine $N_{c}$ for a range of concentrations and ${\rm N_{200}}$. We find that, for $10^{2}<N_{200}<10^{8}$, $C$ varies by a factor \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}2.4 for $c=20$ and \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1.9 for $c=5$. Neglecting this weak dependence, eq. 16 implies that the ratio $r_{\rm conv}/r_{200}$ is approximately independent of mass and that, to first order, $r_{\rm conv}/r_{200}\propto N_{200}^{-1/3}$. As a result, haloes of a given $N_{200}$ will be “converged” roughly to a fixed fraction of their virial radii regardless of mass.

We can use these finding to cast eq. 16 into convenient forms that depend only on particle mass, or mean inter-particle spacing:

[TABLE]

where we have used the fact that $m_{\rm DM}=\Omega_{\rm DM}\,\rho_{\rm crit}(z)\,l(z)^{3}$; $l(z)\equiv l/(1+z)$ is the mean inter-particle spacing in physical units. The latter result, eq. 18, implies that the ratio $r_{\rm conv}(z)/l(z)$ should be largely independent of redshift, halo mass and particle mass; the convergence radius is simply a fixed fraction of the mean inter-particle spacing. We will test these scalings explicitly in later sections.

Figure 1 plots $r_{\rm conv}$ (for $\kappa_{\rm P03}=1$; thick black lines) for NFW profiles with $c=5$ and $c=20$, representative of the vast majority of DM haloes that form in typical cosmological simulations. Due to the weak dependence of $N_{c}$ on $N_{200}$, these curves are only slightly steeper than $r_{\rm conv}/r_{200}\propto N_{200}^{1/3}$. As a result, adopting a fixed softening parameter for cosmological simulations with uniform mass resolution will not necessarily compromise spatial resolution at any mass scale. Take for example the solid blue line in Figure 1, which plots $\epsilon/r_{200}\propto N_{200}^{1/3}$ adopted for the eagle simulation. Here, $\epsilon$ is smaller than $r_{\rm conv}$ by more than a factor of $\sim$2 at all relevant halo masses (but will eventually exceed $r_{\rm conv}$ at very large ${\rm N_{200}}$).

4 Median mass profiles

The most comprehensive attempt to establish the impact of numerical parameters on halo mass profiles has been the work of P03. Their results imply that, provided timestep size, softening and starting redshift are wisely chosen, particle number is the primary factor determining convergence.

One limitation of the work carried out by P03 and N10 is that they were based on simulations of a single $\sim 10^{12}\,{\rm M_{\odot}}$ dark matter halo that was resolved with at least $\sim 10^{4}$ particles, over 300 times that of the poorest resolved systems in typical cosmological runs. Can the conclusions of P03 be extrapolated to haloes with $\sim 10^{2}$ particles, or fewer? Is the P03 radial convergence criterion valid for stacked mass profiles, or for haloes whose masses differ substantially from $\sim 10^{12}\,{\rm M_{\odot}}$? We devote this section to addressing these questions.

4.1 Integration accuracy

The central regions of dark matter haloes are difficult to simulate due to their high densities, which reach many orders of magnitude above the cosmic mean. Crossing times there are $\sim 10^{-3}$ to $10^{-4}\times t_{H}(z)$ implying that thousands of orbits per particle must be integrated to ensure accuracy and reliability. Adding to this, haloes are not smooth and integration errors may accumulate during close encounters between particles. To prevent this, integration must be carried out with a minimum (but a priori unknown) level of precision. In gadget, particles take adaptive timesteps of length $\Delta t=\sqrt{2\eta\epsilon/|\mathbf{a}|}$, where $|\mathbf{a}|$ is the magnitude of the local acceleration, and $\eta$ is a parameter that allows some additional control over the step size. Clearly smaller timesteps are required when $\epsilon$ is small. In the gadget parameter file, $\eta$ is referred to as ${\tt ErrTolIntAccuracy}$ and takes on a default value of 0.025.

Figure 2 plots the average enclosed mass profiles of haloes drawn from our $N_{\rm p}=376^{3}$ run, and highlights the importance of accurate integration. Panels correspond to different mass bins, which were selected so that (from top left to bottom right) $N_{200}\approx 10^{3}$, $10^{4}$, $10^{5}$ and $10^{6}$. The curves show the radii, $r_{\rm N_{p}}$, enclosing a given number of particles (indicated to the right of each curve) and are plotted as a function of the gravitational softening, $\epsilon$. As a guide, the virial radius for each mass bin is indicated by the arrow on the left side of each panel. Connected (blue) circles show results for gadget ’s default integration accuracy parameter, $\tt{ErrTolIntAcc}=0.025$.

For comparison, the vertical green arrows in each panel mark two estimates of minimum softening needed to suppress discreteness effects. The first, $\epsilon_{\rm acc}=r_{200}/\sqrt{{\rm N}_{200}}$, ensures that the maximum stochastic acceleration felt by a particle ($\sim Gm/\epsilon^{2}$) remains smaller than the minimum mean field acceleration of the halo ($\sim GM_{200}/r_{200}^{2}$). The other, $\epsilon_{90}=2\,r_{200}/{\rm N}_{200}$, is less restrictive and is required to prevent large-angle deflections during two-body encounters.

There are a couple points to note in this figure. First, central densities are noticeably suppressed for radii smaller than the “spline” softening length, i.e. r\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\epsilon_{\rm sp}=2.8\times\epsilon (in gadget , pairwise forces become exactly Newtonian for separations larger than $\epsilon_{\rm sp}$). This can be seen by noting that all radii appear to increase slightly once they cross into the grey shaded region, which delineates $r=\epsilon_{\rm sp}$. Spatial resolution is clearly compromised at all radii \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\epsilon_{\rm sp}. More interestingly, however, the central regions of haloes show a clear increase in central density (i.e. the blue lines curve downward) when the softening parameter is reduced below a particular value. The value of $\epsilon$ at which this occurs depends on particle number, with more massive systems exhibiting symptoms at smaller $\epsilon$. Note too that further reducing $\epsilon$ does not result in a denser centre: for very small $\epsilon$, densities are again reduced. The asymptotic central density therefore depends non-monotonically on $\epsilon$, suggesting a numerical origin.

The connected (red) squares show, for a subset of $\epsilon$, the same results but for a series of runs in which $\eta$ was reduced from 0.025 to 0.0025 (this increases the total number of timesteps by a factor of roughly $\sqrt{10}\approx 3.16$). These curves are clearly flatter, suggesting that halo mass profiles are robust to changes in softening across a wide range of masses provided: a) $\epsilon$ is sufficiently small so that r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}\epsilon_{\rm sp}, and b) care is taken to ensure particle orbits are resolved with a sufficient number of timesteps. For the remainder of the paper, all results from runs for which $\epsilon<\epsilon_{\rm fid}$ were carried out with ErrTolIntAcc=0.0025, unless stated otherwise.

These results are qualitatively consistent with those of P03, but differ in the details. These authors report that gadget runs with fixed timestep developed artificially dense “cuspy” centres, where resolution is poor. Based on this, they argue for an $\epsilon$-dependent adaptive timestepping criterion, $\Delta t_{i}=\sqrt{\eta\,\epsilon_{i}/a_{i}}$, where $\eta=0.04$. This is the same criterion adopted for our runs, but with a considerably larger value of $\eta$, perhaps due to the much smaller softening parameters tested in our study. This may indicate the need for a timestepping criterion with a stronger dependence on $\epsilon$. We defer this task to future work.

Provided a sufficient number of timesteps are taken, none of our runs appear to be adversely affected by $\epsilon$ being too small. Indeed, enclosed mass profiles appear stable at all r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}\epsilon even when \epsilon\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\epsilon_{90}, suggesting that hard scattering does not unduly influence the innermost structure of DM haloes, at least for the values of $\epsilon$ tested here. This may be because such scattering events are sufficiently rare (occurring with a cross-section proportional to $\epsilon^{2}$) that their cumulative effects are dynamically unimportant across a Hubble time. Averaged over much longer timescales, scattering-induced cores do develop in idealized simulations of isolated NFW haloes (see, for example, Figure 6 of van den Bosch & Ogiya 2018, who study the secular evolution of an NFW halo for $\approx 60$ Gyr).

These results suggest that the median mass profiles of dark matter haloes are remarkably robust to changes in gravitational softening provided it is not so large that it compromises enclosed masses. This is true even for haloes containing as few as $N_{200}=10^{3}$ particles, and for radii enclosing as few as $50$. This does not, however, imply numerical convergence. Indeed, provided other numerical parameters are carefully chosen, 2-body relaxation places much stronger constraints on convergence (P03), and dense haloes centres – which occupy only a small fraction of their total mass and volume – are highly susceptible to low particle number. The systematic effects of 2-body scattering on the innermost mass profiles of dark matter haloes must therefore be carefully considered when seeking to quantify numerical convergence. As discussed above, a useful tactic for separating converged from unconverged parts of a halo is to identify the radius at which the collisional relaxation time is some multiple $\kappa\equiv t_{\rm rel}/t_{H}(z)$ of a Hubble time (P03). We turn our attention to this in the following subsections.

4.2 Median circular velocity profiles

Figure 3 shows the median circular velocity profiles of haloes in four separate mass bins and at three different resolutions. Grey dashed lines show NFW fits to the profiles of haloes in our $N_{\rm p}=752^{3}$ run (grey lines and circles). Because these curves agree well with the simulated profiles over a large radial range, we can use them to estimate the convergence radius by identifying the point at which the simulated profiles first dip below the theoretical ones by a certain amount. Thick lines (or points in the case of ${\rm N_{p}}=752^{3}$) cover radii for which $V_{c}$ departs from dashed lines by less than 10 per cent; thin lines extend to radii enclosing ${\rm N}\geq 20$ particles.

Note that measured convergence radii are, for a given resolution, only weakly dependent on halo mass (the thick segments or points end at similar radii for a given set of curves). Haloes in our $N_{\rm p}=752^{3}$ run, for example, have convergence radii that differ by at most $\sim$30 per cent across the entire mass range plotted (which corresponds to a factor of $512$ in $N_{200}$ and $8$ in $r_{200}$), and the difference is similar for the other resolutions. Two haloes with the same total number of particles therefore have very different convergence radii depending on their mass. For example, a halo of mass $1.2\times 10^{9}\,{\rm M_{\odot}}$ in our 7523 run has a measured convergence radius of $\approx 0.60\,{\rm kpc}$, whereas $7.5\times 10^{10}\,{\rm M_{\odot}}$ haloes in our 1883 run–resolved with the same number of particles–$r_{\rm conv}\approx 3.5\,{\rm kpc}$, about a factor of 6 larger.

This scaling is indeed expected from eq. 16. Neglecting the weak dependence of concentration on halo mass, two haloes with the same number of particles identified in simulations of different mass resolution, will have convergence radii that resolve comparable fractions of their virial radii, but differ, on average, by a factor $(r_{200}^{\rm A}/r_{200}^{\rm B})\propto(m_{\rm p}^{\rm A}/m_{\rm p}^{\rm B})^{1/3}=({\rm M_{200}^{\rm A}}/{\rm M_{200}^{\rm B}})^{1/3}$, indicating a smaller convergence radius in the higher-resolution run. The downward pointing arrows in Figure 3 mark Power’s convergence radii for each simulation volume and mass bin, which agree well with these empirical estimates. Note that these convergence radii have been approximated using the NFW profiles assuming $\kappa_{\rm P03}=0.177$, smaller than the value advocated by P03. We discuss possible reasons for this difference in Section 6.

4.3 The convergence radius of collisionless cold dark matter haloes

4.3.1 Dependence on halo mass

The convergence radius can also be calculated explicitly by comparing the circular velocity profiles of haloes to those of the same mass but in a higher-resolution simulation. Figures 4 shows, as a function of $N_{200}$, the radius at which the median profiles in our $N_{\rm p}=376^{3}$ (red circles) and 1883 (blue squares) first depart from those in the $N_{\rm p}=752^{3}$ run by more than 10 per cent. Note that only bins containing at least 8 haloes have been included.

For comparison, we also plot the convergence radii expected from eq. 14 as solid lines of corresponding colour (assuming $\kappa_{\rm P03}=0.177$; smaller than the value of $\approx 0.6$ advocated by P03 for $\Delta V_{c}/V_{c}=0.1$), which agree well with the measured values of $r_{\rm conv}$. Our measurements are also recovered well by eq. 15, which is shown using dashed lines in Figure 4 for $\kappa_{\epsilon}=0.177$. Both sets of curves were constructed using NFW haloes that follow the mass-concentration relation of Ludlow et al. (2016), assuming a Planck cosmology, and, for the case of eq. 15, the softening length adopted for each particular run. The shaded region highlights the expected change in $r_{\rm conv}$ brought about by increasing or decreasing $\kappa_{P03}$ by 25 per cent.

Note, however, that eqs. 14 and  15 predict that $r_{\rm conv}$ should depend weakly on mass. While not inconsistent with our numerical results, a much simpler approximation for $r_{\rm conv}$ can be obtained from eq. 18 after we specify a value of $C$. We find $C\approx 2.44$ provides a conservative upper-limit on $r_{\rm conv}$, leading to the following approximation for the P03 convergence radius:

[TABLE]

which is valid for $\Delta V_{c}/V_{c}\approx 0.1$. The dot-dashed horizontal lines in Figure 4 show eq. 20 for our $N_{p}=188^{3}$ (blue) and $376^{3}$ runs (red).

4.3.2 Dependence on gravitational softening

Based on the discussion in Section 3.2, we expect $r_{\rm conv}$ to depend slightly but systematically on $\epsilon$, particularly for haloes resolved with small numbers of particles. Figure 5 plots $r_{\rm conv}$ versus $N_{200}$ for a series of $N_{\rm p}=376^{3}$ runs carried out with different softening lengths. As before, $r_{\rm conv}$ is estimated by comparing the point at which these profiles first depart from those in our highest resolution run ($N_{\rm p}=752^{3}$, $\epsilon=43.75\,{\rm pc}$) by a certain amount. Only bins containing at least 8 haloes are shown. All panels shows results for $\Delta V_{c}/V_{c}=0.1$, with $\epsilon$ decreasing from top to bottom, left to right by a factor of two between consecutive panels. For each value of $\epsilon$ we use filled circles to indicate $N_{200}>2\,r_{200}/r_{\rm conv}$ (or equivalently, r_{\rm conv}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}\epsilon_{90}), which, based on our analytic estimates, should be unaffected by large-angle scattering of particles during close encounters (see Section 3.1). The solid black line in each panel shows the convergence radii expected from eq. 14; dashed lines show eq. 15, which depend explicitly on $\epsilon$ (both assume $\kappa=0.177$). The dot-dashed horizontal lines show eq. 20.

When $\epsilon$ is large, the measured convergence radii scale roughly as $r_{\rm conv}\sim 2\times\epsilon$ (shown as arrows on the right side of each panel), approximately independent of halo mass. Once $\epsilon$ becomes smaller than the analytic estimates of $r_{\rm conv}$ (solid, dashed or dot-dashed lines), the measured values bottom-out and exhibit, at most, a weak dependence on $\epsilon$ and $N_{200}$ thereafter. The weak mass-dependence is described reasonably well by eqs. 14 and  15, but may also be approximated by the much simpler relation, eq. 20, in which $r_{\rm conv}$ is a fixed fraction of the mean inter-particle spacing, regardless of mass. We conclude that, provided $\epsilon$ is negligibly small, eqs. 14 and 15 provide a reasonable upper limit to the values of $r_{\rm conv}$ measured from the median $V_{c}(r)$ profiles of haloes composed of as few as $N_{200}\approx 100$ particles, provided $\kappa_{\rm P03}\approx 0.177$ (for $\Delta V_{c}/V_{c}=0.1$). The dependence of $t_{\rm relax}$ on $\epsilon$ anticipated from eq. 15 adds only a minor correction, but may become increasingly important as $\epsilon$ becomes arbitrarily small.

4.3.3 Dependence on redshift

Convergence radii anticipated from eqs. 14 and 15 depend not only on enclosed particle number, $N(r)$, and gravitational softening, but also on redshift. We show this explicitly in the upper panel of Figure 6, where we plot the convergence radii of haloes identified in one of our $N_{\rm p}=376^{3}$ runs at several different redshifts (this run used $z_{\rm phys}=2.8$ and a maximum physical softening length of $\epsilon_{\rm fid}/2=175\,{\rm pc}$). As in previous Figures, the solid and dashed lines show the analytic estimates of $r_{\rm conv}$ expected from eqs. 14 and 15, which describe the numerical results reasonably well. Convergence radii clearly depend on redshift in a way that is well captured by these simple analytic prescriptions.

Note too the strong redshift-dependence: from $z=0$ to $z=6$, for example, physical convergence radii vary by as much as a factor of $(1+z)=7$ at essentially all mass scales probed by our simulations, consistent with the redshift dependence of $r_{200}(z)$ or $l(z)$. This result may at first seem puzzling, but is indeed expected from eq. 14: haloes that follow a universal NFW mass profile whose concentration depends only weakly on mass and redshift should have convergence radii that scale approximately as $r_{\rm conv}(z)\propto r_{200}(z)/N_{200}^{1/3}$. We show this explicitly in the lower panel of Figure 6, where we have rescaled the convergence radii above by $r_{200}(z)$, and included one of our $N_{\rm p}=188^{3}$ runs (blue points). All curves now follow a similar scaling, implying that the co-moving convergence radius is largely independent of redshift. The dot-dashed black line in the lower panels shows eq. 20, for which $r_{\rm conv}$ is a fixed fraction of the mean inter-particle spacing: $r_{\rm conv}/l=0.055$. This simple approximation describes our numerical results well, but can be improved slightly using $r_{\rm conv}/r_{200}=1.1\times N_{200}^{-0.4}$ (heavy solid line in the lower panels), which is slightly steeper than $N^{-1/3}$. These convergence radii are similar to, but slightly larger than, the “‘optimal” softening for NFW haloes advocated by van den Bosch & Ogiya (2018). For $N_{200}$ spanning $\sim 10^{2}$ to $\sim 10^{7}$, $\epsilon_{\rm opt}$ is a factor of 2 to 3 smaller than these estimates of $r_{\rm conv}$, and should therefore not compromise spatial resolution. In addition, since $\epsilon_{\rm opt}\propto N_{200}^{-1/3}$, the ratio $\epsilon_{\rm opt}/m_{\rm DM}^{1/3}$ remains fixed for all $N_{200}$: the softening is optimal at all masses, regardless of $N_{200}$, making it a potentially desirable choice for cosmological simulations of equal mass particles.

4.3.4 Dependence on $\kappa$

In previous sections we estimated convergence radii in our simulations by comparing the median circular velocity profiles of haloes in our $N_{\rm p}=376^{3}$ simulation with those of the same virial mass but in a higher-resolution simulation ($N_{\rm p}=752^{3}$). The radius within which profiles deviate by more than 10 per cent marked $r_{\rm conv}$. For analytic estimates of $r_{\rm conv}$ based on eqs. 14 and 15 this corresponds to particular values of $\kappa$, about $0.177$. However, as discussed in detail in N10, better convergence can be obtained for higher values of $\kappa$, which occurs at larger radii where relaxation times are substantially longer.

In Figure 7 we plot three separate estimates of convergence radii, corresponding to fractional departures between median $V_{c}(r)$ profiles in our low and high-resolution runs of $\Delta V_{c}/V_{c}=0.03$, 0.1 and 0.15. Better convergence is obtained at larger radii, and requires larger values of $\kappa$: the solid and dashed lines show convergence radii estimated from eqs. 14 and 15, respectively; blue, orange and green denote $\kappa=0.566$, 0.177 and 0.106, respectively (as before, the shaded region indicates $\pm 25$ per cent in $\kappa_{\rm P03}$). Horizontal dot-dashed lines indicate fixed fractions of the mean inter-particle distance, corresponding to $r_{\rm conv}/l=0.040$, 0.055 and 0.10 for $\Delta V_{c}/V_{c}=0.03$, 0.1 and 0.15, respectively.

Figure 8 summarizes a number of previous results. Here we plot, in the top panels, the residual difference in the circular velocity profiles between haloes in our low- and highest-resolution runs as a function of $\kappa(r)=t_{\rm rel}(r)/t_{H}(z)$. The left hand panel shows results for $N_{\rm p}=376^{3}$ (red) and 1883 (blue); both runs used a ($z=0$) softening length of $\epsilon=43.75\,{\rm pc}$. The panel on the right compares several $N_{\rm p}=376^{3}$ runs using different softening. In all cases, residuals are calculated with respect to our $N_{\rm p}=752^{3}$ ($\epsilon=43.75\,{\rm pc}$) in 40 equally-spaced bins of $\log\,M_{200}$ spanning the range $100\,m_{\rm DM}$ to $10^{12}\,{\rm M}_{\odot}$. White circles and squares mark the empirical results of P03 and N10, respectively; both are more conservative than ours (which supports the conclusions of Zhang et al., 2018).

As expected from previous plots, deviations in $V_{c}$ are largely independent of both mass and force resolution, but correlate strongly with the enclosed relaxation timescale $\kappa(r)$. Overall, the convergence of the median circular velocity profiles can be approximated reasonably well by

[TABLE]

where we have defined $\psi\equiv\log\kappa+d$, and $(a,b,c,d)=(-0.4,-0.6,-0.55,0.95)$ (heavy black line in the upper panels of Figure 8).

In the lower panels of Figure 8 we plot $\Delta{\rm V_{c}/V_{c}}$ for the same runs, but now as a function of halo-centric distance normalized by the mean inter-particle separation, $l$. Convergence in circular velocity is achieved at spatial scales that roughly correspond to fixed fractions of $l$, which provides a much simpler convergence criterion than that advocated by P03. A conservative upper limit to the convergence radius as a function of $\Delta{\rm V_{c}/V_{c}}$ can again be approximated by eq. 21, but with modified parameters: $\psi\equiv\log(r/l)+d$ and $(a,b,c,d)=(-1.96,-0.76,-0.52,1.49)$ (thick black line in the lower panels). The outsized color points in Figure 8 correspond to the curves in Figure 7.

4.4 Scaling relations

These results help clarify why structural scaling relations for dark matter haloes converge even for systems resolved with only a few hundred particles, a result we show explicitly in Figure 9. Here we plot, as a function of $M_{200}$, the ratio of the virial radius to the radius $r_{X\%}$ enclosing X per cent of the halo mass (the curve labelled “50%”, for example, is the half-mass radius–virial mass relation), for all haloes resolved with $N_{200}\geq 32$ particles. As in previous figures, different colours correspond to different resolutions. Points connected by thick lines in Figure 9 correspond to runs carried out with fixed $\epsilon=43.75\,{\rm pc}$, and the thin lines to the “fiducial” softening, which is a factor of 4 to 16 times larger, depending on resolution. The faint diagonal lines show the P03 convergence radius for $\kappa_{\rm P03}=0.177$ (solid) and assuming $r_{\rm conv}/l=0.055$ (dotted). These provide a good indication of the mass scale above which the scaling relations are converged. Note too that, as expected, the converged relations are largely independent of $\epsilon$. Dashed lines show the expected trends assuming NFW profiles and the $c(M)$ relation of Ludlow et al. (2016). For comparison we also plot the expected concentration, $c=r_{200}/r_{-2}$, and the ratio $r_{200}/r_{\rm max}$ using the same model, as well as radii that enclose fixed overdensities of $\Delta=500$ and 2500.

Resolving the innermost structure of haloes is clearly challenging. Resolving $r_{-2}$ with a systematic bias in ${\rm V_{c}}(r_{-2})$ of less than $10$ per cent (assuming eq. 14 with $\kappa_{\rm P03}=0.177$), for example, requires $N_{200}\approx 844,\,1083$, and $1353$ particles for our $N_{\rm p}=188^{3}$, $376^{3}$ and $752^{3}$ runs, respectively. Resolving $r_{\rm max}\approx 2.2\times r_{-2}$ is less demanding, requiring only $N_{200}\approx 170,\,214,$ and $266$ for the same respective $N_{\rm p}$.

We show this explicitly in the left-hand panel of Figure 10, where we plot the $V_{\rm max}-r_{\rm max}$ relations for haloes in runs of different resolution. Thick curves correspond to the median $r_{\rm max}$ and $V_{\rm max}$ in bins of $M_{200}$, and extend down to a lower mass limit corresponding to $N_{200}=100$. We adopt the same colour scheme as in previous figures and use solid lines for our fiducial runs, and dashed lines for runs where $\epsilon$ was kept fixed at $43.75\,{\rm pc}$. Note that curves corresponding to a given mass resolution begin to converge when $N_{200}$ exceeds the lower limits provided above. Convergence is not perfect, however, as systematic differences in $V_{\rm max}$ and $r_{\rm max}$ of order 10 per cent are expected at these mass scales. Solid (faint) lines of corresponding colour, for example, show the $r_{\rm max}-V_{\rm max}$ relation for haloes whose mass is kept fixed at those lower limits (the relevant curves are labelled $r_{\rm max}-r_{\rm conv}$). Dotted lines show the analogous relations for convergence radii equal to $r_{\rm conv}=0.055\,l$.

Fixing $\epsilon$ at $43.75\,{\rm pc}$, well below the fiducial value, appears to improve convergence between runs of different mass resolution even slightly *below * these mass scales. Indeed, at first glance, convergence even seems better than expected from limits imposed by 2-body relaxation. This result, however, is fortuitous, as shown in the right-hand of Figure 10 where $r_{\rm max}-V_{\rm max}$ relations are plotted for our $N_{\rm p}=376^{3}$ runs for a variety of $\epsilon$. It is clear from this figure that a r_{\rm max}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}r_{\rm conv} is a requirement for convergence in the median relations (i.e. convergence is only achieved to the right of the solid or dotted grey lines labelled $r_{\rm max}=r_{\rm conv}$).

4.5 Convergence Requirements for cosmological simulations

The results of previous sections suggest that we can use eq. 14 or  15 to place restrictions on the minimum number of particles required to reach a desired spatial resolution in halo centres. The left panel of Figure 11 plots the total number of particles, $N_{200}$, required to resolve the innermost $r_{\rm conv}=1\,{\rm kpc}$ (green), 0.1 ${\rm kpc}$ (orange) and 0.01 ${\rm kpc}$ (blue) of an NFW halo as a function of its virial mass $M_{200}$. Solid and dashed lines correspond to $\kappa_{\rm P03}=0.177$ and 0.566 (using eq. 14), respectively. All curves assume the mass-concentration relation predicted by the model of Ludlow et al. (2016). This plot highlights the difficulty of resolving the inner regions of dark matter haloes. For example, to resolve the central $1\,{\rm kpc}$ of a $10^{15}\,{\rm M_{\odot}}$ galaxy cluster with $10$ per cent accuracy (corresponding to $\approx 0.06$ per cent of $r_{200}$) requires sampling its virial mass with $>10^{8.4}$ particles, comparable to the particle number required to resolve the innermost $100\,{\rm pc}$ of Milky Way mass haloes. Note too that achieving a spatial resolution of $\sim 10\,{\rm pc}$ for a simulated halo of mass comparable to that of the Milky Way would require $N_{200}\sim 10^{11}$ particles, far higher than any cosmological simulation published to date. Note, however, that these scales are baryon-dominated in hydrodynamical simulations, and the relevance of these criteria are not obvious in that case.

The middle panel of Figure 11 shows the convergence radius (computed using eq. 14 assuming an NFW profile and $\kappa_{\rm P03}=0.177$) as a function of halo mass expected for simulations of different uniform mass resolution (burgundy lines), for which $m_{\rm DM}$ varies by successive factors of 8. Heavy lines highlight several particular values of $m_{\rm DM}$. Note that for $m=1.2\times 10^{7}\,M_{\odot}$ (comparable to the particle mass in the $100\,{\rm Mpc}$, $N_{\rm p}=1504^{3}$ eagle simulation), convergence radii vary from $\approx 3\,{\rm kpc}$ for $10^{10}\,{\rm M}_{\odot}$ haloes, to $\approx 2\,{\rm kpc}$ at ${\rm M_{200}}=10^{15}\,{\rm M}_{\odot}$, equivalent to roughly 3 to 4 fiducial softening lengths. Targeting convergence radii below \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100\,{\rm pc} for haloes with virial masses M_{200}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{8}\,M_{\odot} requires dark matter particle masses of only a few thousand ${\rm M}_{\odot}$; achieving \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\,{\rm pc} resolution requires $m_{\rm DM}\approx{\rm M}_{\odot}$.

Dashed blue lines show, for comparison, the convergence radii as a function of mass for fixed values of $N_{200}$, ranging from 20 particles (thick dashed line) to $N_{200}=10^{9}$. This is comparable to the highest-resolution dark matter-only simulations carried out to date – the Aquarius (Springel et al., 2008) and Ghalo (Stadel et al., 2008) simulations – which achieve a maximum spatial resolution of order $100\,{\rm pc}$ in a Milky Way mass dark matter halo. These results imply that the current state-of-the-art in cosmological hydrodynamical simulations (e.g. eagle, Illustris, Apostle, FIRE, Auriga) are not fully converged on scales relevant for galaxy formation.

The right-most panel of Figure 11 plots the number of particles within $r_{\rm conv}$ as a function of $N_{200}$ for several different values of the (NFW) concentration (we assumed $\kappa_{\rm P03}=0.177$ and chose values of $c$ that span the extremes of concentrations measured in typical cosmological simulations). Once a mass profile and concentration have been chosen, $N(<r_{\rm conv})$ depends only on $N_{200}$. Haloes with $N_{200}\approx 10^{3}$ will typically have between $\approx 30$ and 100 particles within $r_{\rm conv}$ for $c=3$ and 20, respectively; haloes with $N_{200}\approx 10^{6}$ will have between $\approx 260$ and 1500, respectively.

4.6 Section summary

Before moving on, we first summarize the results of this section:

• •

Provided a sufficient number of timesteps are taken, the median mass profiles of simulated haloes are independent of softening at virtually all radii similar to or greater than the spline softening length, r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}\epsilon_{\rm sp}. This confirms the findings of P03, but extends their results to a broader range of halo mass, and to more extreme values of $\epsilon$. This holds true even for runs with very small softening and for haloes resolved with only $\sim 10^{3}$ particles, despite the fact that strong discreteness effects were naively expected for \epsilon\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}r_{200}/\sqrt{N_{200}}. Indeed, the enclosed mass profiles of dark matter haloes are remarkably robust to changes in $\epsilon$, even for much smaller values of $\epsilon$.

• •

Nevertheless, the fact that the median mass profiles of haloes are largely insensitive to $\epsilon$ says little about the radial range over which they can be considered reliably resolved. For a given DM halo the minimum resolved radius, $r_{\rm conv}$, depends primarily on total particle number, $N_{200}$, and roughly coincides with the radius within which the enclosed 2-body relaxation time first exceeds the Hubble time by some factor $\kappa$. The precise value of $\kappa$ dictates the level of convergence: we find that for $\kappa\approx 0.177$, median circular velocity profiles have converged to within 10 per cent; 3 per cent convergence in $V_{c}(r)$ requires $\kappa\approx 0.566$. This is true regardless of $m_{\rm DM}$ and of $\epsilon$, provided the condition \epsilon\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}r_{\rm conv} is met. Note that these values of $\kappa$ are smaller than those advocated by N10 and P03 for similar levels of convergence (e.g. $\kappa=7$ for 2.5 per cent convergence, and $\kappa=0.6$ for 10 per cent), suggesting that median halo mass profiles may converge more easily than those of individual systems.

• •

The convergence radius scales roughly as $r_{\rm conv}/r_{200}\propto N_{200}^{-1/3}$ (see eq. 16), implying $r_{\rm conv}/l\approx{\rm constant}$, independent of redshift, halo mass and particle mass. Indeed, for $\Delta V_{c}/V_{c}=0.1$ ($\kappa_{\rm P03}=0.177$) we find that $r_{\rm conv}(z)=0.055\times l(z)$ describes our numerical data about as well as the P03 criterion, typically to within 20 per cent. A better approximation can be obtained with a slightly steeper dependence on $N_{200}$: $r_{\rm conv}/r_{200}=1.1\times N_{200}^{-0.4}$ (see Figure 6).

• •

Cosmological simulations should adopt softening lengths at least a factor of 2 smaller than $r_{\rm conv}$ in order to ensure that biased force estimates do not compromise their spatial resolution. Recently, van den Bosch & Ogiya (2018) suggested that the “optimal” softening for NFW haloes scales with particle number approximately as $\epsilon_{\rm opt}/r_{200}=0.005\times(N_{200}/10^{5})^{-1/3}$, which is a factor of 2 to 4 smaller than $r_{\rm conv}$ over virtually all values of $N_{200}$ relevant for cosmological simulations ($10^{2}\leq N_{200}\leq 10^{8}$), and is “optimal” regardless of halo mass. Indeed, their results imply that $\epsilon_{\rm opt}/l\approx 0.017$, comparable to values adopted by the majority of recent large-scale cosmological simulations.

5 Halo mass functions

Accurately characterizing the dark matter halo mass function, $n(M)$, either from simulation or theory, is necessary for several reasons. At high mass, for example, the shape and time evolution of $n(M)$ encodes clues that provide important constraints on cosmological models (e.g. Eke et al., 1996). The mass function is also an important probe of dark matter, since many particle candidates predict strong, scale-dependent deviations from the expectations of the canonical cold dark matter model (see, e.g. Schneider, 2015; Angulo et al., 2013). Accurate predictions for halo mass functions are also required for galaxy formation theory, since galaxies are assumed to form in halo centres via dissipative collapse of primordial and recycled gas (see, e.g., White & Frenk, 1991, and a number of subsequent works). The abundance of dark matter haloes is therefore closely related to the observed luminosity function of galaxies, and is a central aspect of theoretical models that attempt to reproduce observed galaxy number densities.

The latter point is particularly important for hydrodynamical simulations that possess star-forming haloes ($\;\buildrel>\over{\sim}\;$ a few $\times 10^{9}\,{\rm M_{\odot}}$) close their resolution limit. eagle is an example of a simulation in which haloes of mass $\sim 10^{9}\,{\rm M}_{\odot}$ are resolved with only a few hundred dark matter particles. The abundance and structure of these haloes, and consequently the statistics of the first generation of low-mass galaxies, are likely subject to numerical artifact. In this section we quantify the robustness of the halo mass function to changes in particle mass and in gravitational softening length, focusing particularly on the lowest mass haloes.

5.1 Dependence on halo mass definition

Before studying the sensitivity of halo mass functions to numerical parameters, it is useful to examine their dependence on how mass is defined, in order to have a useful gauge for later on. Figure 12 shows the cumulative mass functions of haloes in our $N_{\rm p}=752^{3}$ ($\epsilon=43.75\,\rm{pc}$) run for several possibilities. The connected points correspond to FoF masses, solid blue and green lines to the overdensity masses ${\rm M}_{200}$ and ${\rm M}_{vir}$, respectively. Red and yellow curves use the bound halo mass either excluding (red) or including (yellow) the contribution of its substructure population (see Section 2.2 for details). Residuals in the lower panel are computed with respect to the FoF mass definition.

All mass functions have a similar shape, which can be approximated by a power-law over the mass range covered by the run. The residuals, however, betray more clearly the differences. Based on mass definition alone, $n(M)$ may differ by of order a few to $\approx$20 per cent for haloes resolved with \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}100 particles, with differences becoming even larger towards lower mass. FoF masses generally result in the highest number densities at fixed $\rm M_{\rm FoF}$, and those based on $\rm M_{200}$ the lowest. Note that even relatively well-resolved haloes have number densities that differ by roughly 10 per cent between $\rm M_{\rm FoF}$ and $\rm M_{200}$ (corresponding to an average systematic mass difference of roughly 15 to 20 per cent).

We next consider the sensitivity of halo mass functions to numerical parameters, focusing on mass and force resolution. As we will see, systematic differences brought about by varying numerical parameters are virtually always smaller than those associated with mass definition.

5.2 Dependence on particle number and softening

The upper panels of Figure 13 show the cumulative mass functions for ${\rm M_{FoF}}$ (left) and ${\rm M_{200}}$ (right) in our $N_{\rm p}=752^{3}$, 3763 and 1883 boxes. Solid lines correspond to runs carried out with the fiducial softening parameter, $\epsilon/b\approx 0.011$ (at $z=0$); dashed lines are used for runs in which $\epsilon$ was kept fixed at $43.75\,{\rm pc}$ (physical, also at $z=0$), regardless of $m_{\rm DM}$. The upper residual panel plots the departure of each fiducial run from the black dashed line (i.e. from the run with the highest mass resolution and smallest force softening). Differences are small, typically \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5 per cent for haloes containing {\rm N}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}100 particles (indicated using downward-pointing arrows), and \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10 per cent at all masses resolved by our simulation. The lower residual panels compare results for fixed-softening runs (we use $\epsilon=43.75\,{\rm pc}$, corresponding to $\epsilon/\epsilon_{\rm fid}=1/4$, $1/8$ and $1/16$ for $N_{\rm p}=752^{3}$, 3763 and 1883, respectively). Albeit slightly, this choice of softening improves agreement at low FoF masses, but has the opposite effect for $\rm M_{200}$. Overall, the $\rm{M_{FoF}}$ mass functions appear more stable to changes in $N_{\rm p}$ than those using ${\rm M_{200}}$. Nevertheless, it is worth emphasizing that, in both cases, departures remain small compared to variations arising from different mass definitions.

The lower two panels of Figure 13 show, for $N_{\rm p}=376^{3}$, the $\rm M_{\rm FoF}$ (left) and $\rm M_{200}$ (right) mass functions for a large range of $\epsilon$ (5.46\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\epsilon/[{\rm pc}]\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5.6\times 10^{3}). The dashed black lines are the same as in the upper panels, and correspond to the run with $N_{\rm p}=752^{3}$, $\epsilon=43.75\,{\rm pc}$; residuals in the lower panels are computed with respect to these curves. Provided \epsilon/\epsilon_{\rm fid}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}8, differences due to softening are most pronounced for the lowest mass haloes. Indeed, for a broad range of softening, 4\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}\epsilon/\epsilon_{\rm fid}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}1/64, both sets of mass functions agree with that of the higher-resolution run to within about $10$ per cent for haloes resolved with more than $\sim$100 particles, although deviations depend systematically (and non-monotonically; see Figure 14 below) on $\epsilon$.

Agreement between the halo mass functions depends on mass and force resolution in different ways depending on how masses are defined. As discussed by Tinker et al. (2008) (see, also, Lukić et al. 2009), SO masses, such as $\rm M_{200}$, are sensitive to the integrated mass profile within a fixed overdensity, whereas FoF masses are measured within iso-density contours, and are less sensitive to the precise distribution of mass within them. We therefore expect FoF mass functions to be more robust to changes in numerical parameters than those based on spherical overdensities, a claim that is supported by the results presented in Figure 13. Indeed, Tinker et al. (2008) show that convergence in SO mass functions worsens considerably with increasing overdensity.

Lukić et al. (2007) and Power et al. 2016 showed that, unsurprisingly, halo mass functions are heavily suppressed on scales below which $\epsilon$ exceeds the halo virial radius; none of our runs are in that regime (the largest softening we test, $\epsilon=5.6\,{\rm kpc}$, is equal to the virial radius of a halo of mass ${\rm M_{200}}\approx 2\times 10^{7}\,{\rm M_{\odot}}$, well below the mass limit of our halo finder). Our runs address a more subtle question: what impact does softening have on the abundance of haloes when $\epsilon\ll r_{200}$?

Figure 14 shows the softening dependence of the total number density of haloes in our $N_{p}=376^{3}$ runs based on $\rm M_{\rm FoF}$ and $\rm M_{\rm 200}$ masses. From top to bottom, the different panels show results for $N\geq 32$, 100, 256 and 800 particles, respectively; the corresponding mass thresholds are indicated in each panel. The filled and open squares correspond to $\rm M_{\rm FoF}$ and $\rm M_{\rm 200}$, respectively; symbols are colour-coded by $\epsilon$ as in Figure 13. For comparison, filled and open circles show the number densities of haloes above the same mass thresholds in the $N_{\rm p}=752^{3}$ runs. The grey shaded regions centred on these circles indicate a 5 per cent change in total abundance. When including the lowest-mass haloes ($N\leq 100$), the total number density grows gradually until $\epsilon\approx 4\times\epsilon_{\rm fid}$, where it peaks, before decreasing rapidly for larger $\epsilon$. For eagle’s fiducial softening length halo abundances are within a few per cent of the maximum attained for any $\epsilon$. However, for both mass definitions, the total abundance of haloes across all $\epsilon$ varies by \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}5 per cent for $N\leq 100$, suggesting that per cent-level convergence in halo mass functions demands higher particle numbers.

What is the true abundance of haloes above each mass threshold? It is tempting to assume that convergence is achieved when the results of low- and high-resolution runs agree. If this were the case, then convergence in the abundance of haloes containing a few dozen particles is clearly ruled out (in the upper panel of Figure 14 these runs never overlap). It is important to note, however, that the $\epsilon-$dependence of $n(M)$ weakens as $N$ increases, and agreement between runs of different mass resolution improves. For $N\geq 256$ (second panel from the bottom), a weak dependence on $\epsilon$ is still evident at the level of $\approx 5$ per cent in our $N_{\rm p}=376^{3}$ run, but has largely disappeared for $N\geq 800$ (bottom panel), where the total abundances of haloes in both runs agree (within the Poisson noise) for all values of $\epsilon$ studied (except, perhaps, for the largest value and $N_{\rm p}=376^{3}$).

If $N\geq 800$ is sufficient to achieve convergence, then results from our $N_{\rm p}=752^{3}$ runs represent robust measurements of $n(M)$ in all but the upper-most panel (haloes of fixed mass are resolved with 8 times the number of particles in this run, so have $N\geq 800$ in each of the lower three panels). In all cases, our low-resolution runs are able to reproduce these abundances, but only for particular values of $\epsilon$, which differ depending on how mass is defined (with FoF masses preferring small softening, \epsilon\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\epsilon_{\rm fid}/8, and $\rm M_{200}$ preferring large values, \epsilon\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}2-4\,\epsilon_{\rm fid}). We conclude that, regardless of mass definition, halo mass functions exhibit a small ($\approx 5$ per cent) but measurable dependence on $\epsilon$ for haloes resolved with \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}800 particles.

5.3 Dependence on redshift evolution of $\epsilon$

The gravitational softening length, initially fixed in comoving coordinates, reaches a maximum physical value at $z_{\rm phys}=2.8$, after which it remains constant in proper coordinates. Other cosmological simulations often use a fixed comoving softening at all $z$, whereas some opt for fixed proper softening lengths. What effect, if any, does this have on the halo mass function? In order to find out, two additional “fiducial” ($N_{\rm p}=376^{3}$) runs were carried out using $z_{\rm phys}=10$ and 0 (i.e., fixed in comoving units at all times), each reaching a $z=0$ softening parameter of $\epsilon_{\rm fid}=350\,{\rm pc}$. The upper panel of Figure 15 shows the redshift evolution of the total number of haloes above $32$ and $100-$particle thresholds. As above, masses are computed at all $z$ using both $M_{\rm FoF}$ (solid lines) and $M_{\rm 200}$ (dashed) definitions. The connected (blue) circles in this plot show the results for $z_{\rm phys}=2.8$; other curves correspond to $z_{\rm phys}=10$ (green diamonds) and 0 (red squares).

Overall, the results are similar, with the largest differences being limited to the lowest-mass haloes. The middle and bottom panels of Figure 15 show the residuals with respect to the run with $z_{\rm phys}=0$ (for clarity, the comparison is limited to ${\rm M_{200}}$ masses in this plot). For $N_{200}\geq 100$ (lower panel), differences are never larger than $\approx 2$ per cent (highlighted by the grey shaded region). For $N_{200}\geq 32$ differences are still small, \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5 per cent, but systematic: higher $z_{\rm phys}$ corresponds to higher $n(M)$. This is qualitatively consistent with the results in Figure 14, which indicates that increasing $\epsilon$ results in a slight boost in the numbers of low-mass haloes, at least up to a point. For a given $z=0$ softening length, higher $z_{\rm phys}$ implies larger physical softening lengths at $z>z_{\rm phys}$, and also enhances slightly the numbers of low-mass haloes. Note as well the slight boost in $n(M)$ (of order a couple per cent) after $z_{\rm phys}$. These subtle changes in the abundance of low-mass haloes may impact the star formation histories of low-mass galaxies that inhabit them, provided they are sufficiently massive to promote efficient gas cooling.

5.4 Section summary

The results of this section can be summarized as follows:

• •

Dark matter halo mass functions may differ by as much as 20 per cent for haloes resolved with $\geq$100 particles depending on how mass is defined. At fixed mass, $M_{200}$ masses yield the lowest overall number densities and $M_{\rm FoF}$ the highest. Comparing haloes at fixed number density, this suggests that, on average, FoF masses exceed those based on $M_{200}$ by a factor of $\approx 1.15-1.2$, though a more detailed comparison is required to properly assess the systematics (see, e.g., Tinker et al., 2008, for a thorough discussion).

• •

For haloes containing $\;\buildrel>\over{\sim}\;$100 particles, FoF mass functions converge to better than $\approx$5 per cent for runs carried out with our “fiducial” softening, $\epsilon_{\rm fid}/l=0.011$ at $z=0$, regardless of particle number $N_{\rm p}$ (note $\epsilon_{\rm fid}=700$, 350, and $175\,{\rm pc}$ for $N_{\rm p}=188^{3}$, 3763 and 7523, respectively), and to $\;\buildrel<\over{\sim}\;$3 per cent when $\epsilon=43.75\,{\rm pc}$ was kept fixed for all $N_{\rm p}$. For haloes resolved with $\geq 32$ particles, corresponding to the resolution limit of subfind, FoF mass functions converge to within $\approx 10$ per cent for fiducial softening values, and to within $\approx 5$ per cent for $\epsilon=43.75\,{\rm pc}$. The overall trend is such that lower mass resolution results in a systematic increase in the total number of FoF haloes, driven mainly by a slight increase in the number of poorly resolved systems. For $M_{200}$ the converse is true: runs of lower resolution produce systematically fewer haloes. Differences are small but systematic, reaching $\approx 5$ per cent in our fiducial runs for haloes resolved with at least $N_{200}\geq 100$, and $\approx 10$ per cent for runs with $\epsilon=43.75\,{\rm pc}$.

• •

At fixed mass resolution (corresponding to $N_{\rm p}=376^{3}$) $M_{\rm FoF}$ mass functions converge to within $\approx 6$ (15) per cent for $N_{\rm FoF}\geq 100$ (32), for a range of $\epsilon$ spanning nearly a factor of $10^{3}$ ($5.46\,{\rm pc}$ to $5.6\,{\rm kpc}$). For $M_{200}$, we obtain convergence at the 11 and 15 per cent level for N_{\rm 200}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}100 and \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}32, respectively, although larger differences are found for the most extreme values of $\epsilon$ tested. There is also a systematic trend: increasing $\epsilon$ results in a gradual increase in the abundance of haloes composed of \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100 particles, provided it remains smaller than about $8\times\epsilon_{\rm fid}\approx 0.09\,l$, where the total abundance of haloes peaks.

• •

Changing the redshift $z_{\rm phys}$ below which the softening parameter remains fixed in proper units (rather than comoving) can have a comparable affect on the abundance of low-mass \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100-particle haloes. When $\epsilon$ is fixed at $z=0$ a higher $z_{\rm phys}$ implies larger physical softening lengths for all $z>z_{\rm phys}$. As we have seen above, larger softening lengths tend to enhance slightly the formation of low-$N$ haloes, provided $\epsilon$ does not become too large. For example, for $N_{\rm p}=376^{3}$ and $\epsilon(z=0)=350\,{\rm pc}$, increasing $z_{\rm phys}$ from 0 to 10 results in a $\;\buildrel<\over{\sim}\;$2 per cent increase in the abundance of haloes with $N_{200}\geq 100$; for $N_{200}\geq 32$ deviations are not larger than $\approx 5$ per cent.

Overall, our results suggest that halo mass functions are a robust result of N-body simulations once a definition of halo mass has been specified. Variations in halo abundances with numerical parameters tend to be restricted to poorly resolved systems containing fewer than $100$ particles. Yet, however slight, these changes may have a noticeable impact on the first generation of star formation if these haloes happen to have masses comparable to the threshold for efficient gas cooling. Calibrating sub-grid models at fixed $N_{\rm p}$ and $\epsilon$ may help mitigate any non-physical effects brought about by star formation in this first generation of poorly resolved haloes, but may not easily adapt to increasing or decreasing resolution. As a result, simulations that adopt numerical parameters that differ from those for which the models were calibrated may yield noticeably different star formation histories and/or galaxy properties. We will address these issues in a companion paper.

6 Conclusions

We have carried out a systematic convergence study of the median and statistical properties of DM haloes in fully cosmological, dark matter-only simulations. Unlike previous work, which targeted single haloes, we focused our analysis on their median spherically-averaged density profiles as a function of mass and on several mass-dependent structural scaling relations. After verifying the need for fine timesteps to resolve halo centres, we tested the sensitivity of halo profiles to total particle number, $N_{\rm p}$, and force softening length, $\epsilon$. We also revisited the calculation of the convergence radius originally provided by P03, and derived its explicit (but weak) dependence on the gravitational softening (eq. 15).

In addition to mass profiles, we also studied convergence in the abundance of haloes as a function of mass, focusing on runs carried out with different particle numbers and redshift-dependent softening lengths. Our main results can be summarized as follows:

Softening and 2-body relaxation: Softening does not significantly affect 2-body relaxation times, which are primarily driven by particle number. This is contrary to common belief, despite being well documented in the literature (e.g. Huang et al., 1993; Theis, 1998; Dehnen, 2001). The result may at first seem puzzling because encounters between particles give rise to velocity perturbations that scale approximately as $\delta v\sim G\,m_{\rm DM}/(b\,v)$, increasing significantly for encounters with small impact parameter, $b$. Nevertheless, close encounters are rare and, as pointed out by Chandrasekhar (1942) and Spitzer & Hart (1971), the large numbers of distant encounters dominate the cumulative effect of perturbations: relaxation is driven by discreteness on large scales, where most of the mass is (see, e.g., Hernquist & Barnes, 1990). Softening does, however, suppress small-scale “collisions” between particles: for Plummer softening interactions between particles become increasingly unimportant for separations \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\epsilon/\sqrt{2}. Softening therefore serves the purpose of smoothing the matter distribution on small scales, allowing force estimates to more faithfully represent a continuous matter field. 2. 2.

Convergence of the median mass profiles of CDM haloes: Provided $\epsilon$ and timestep size are appropriately chosen, 2-body relaxation imposes a strict and well-defined lower limit on the spatial resolution of collisionless CDM haloes (Figure 4). Convergence in mass profiles, for example, can only be achieved at radii beyond which 2-body relaxation times are a sizeable fraction of a Hubble time, or longer. As a result, convergence radii are primarily determined by the enclosed particle number, though analytic arguments suggest a weak (approximately logarithmic) dependence on softening. Softening appears to compromise the spatial resolution of simulations only if it is larger than the convergence radius dictated by 2-body scattering: in this case, our results suggest that $r_{\rm conv}\approx 2\times\epsilon$ (Figure 5).

Particle interactions tend to drive the local velocity distribution towards a Maxwellian and, occasionally, impart velocities exceeding the local escape speed. The net effect of 2-body relaxation is therefore to monotonically suppress the central densities of haloes, allowing convergence radii to be determined empirically by comparing runs of different (mass) resolution. A simple but accurate description of our measured convergence radii can be obtained using eq. 14 or 15, where $\kappa\equiv t_{\rm relax}/t_{\rm H}$; better convergence is achieved for higher values of $\kappa$. Circular velocities, for example, converge to within $\approx 3$, 10 and 20 per cent for $\kappa=0.566$, 0.177 and 0.106, respectively (Figure 7; values of $\kappa$ for arbitrary levels of convergence can be approximated using eq. 21).

These results, valid at all redshifts (Figure 6), are qualitatively consistent with those of P03 and N10, but differ in the details. P03, for example, find that $10$ per cent convergence in $V_{c}(r)$ requires $\kappa\approx 0.6$, while N10 find $\kappa=7.5$ yields $\Delta V_{c}/V_{c}\approx 0.025$; both are larger than the values of $\kappa$ we advocate for the same level of convergence ($\approx 0.2$ and $\approx 0.6$, respectively). It is worth emphasizing, however, that both previous studies focused on single dark matter haloes of fixed virial mass, whereas our results apply to median mass profiles that span a broad range of halo masses. Differences between our results and theirs may arise as a result of sampling haloes possessing a wide range of concentrations at fixed mass: since $r_{\rm conv}$ (and its relation to $\kappa$) depends explicitly on profile shape, it will necessarily vary from halo to halo. Another possibility is that, by averaging, we smooth-out idiosyncrasies of individual systems–such as asymmetric shapes, or locally dominant substructure–thereby minimizing the effects of collective relaxation. The convergence criteria advocated by P03 and N10 are therefore more conservative than ours, a result echoed by the recent work of Zhang et al. (2018).

At fixed mass resolution, the convergence radii anticipated by eq. 14 exhibit a weak dependence on $N_{200}$, varying by about a factor of two for NFW haloes with $N_{200}$ ranging from $\sim 10^{2}$ to $\sim 10^{8}$. Indeed, the measured median convergence radii of haloes in our simulations are compatible with a much simpler approximation in which $r_{\rm conv}$ is simply a fixed fraction of the mean inter-particle spacing, $l$. For example, $<10$ per cent convergence in circular velocity is obtained approximately at a radius $r_{\rm conv}(z)=0.055\times l(z)$ (eq. 20; Figure 4). Softening lengths for cosmological simulations should be chosen with this in mind. 3. 3.

The optimal softening for cosmological simulations: We note that the “optimal” softening length, $\epsilon_{\rm opt}/r_{200}=0.005\times(N_{200}/10^{5})^{-1/3}$, advocated by van den Bosch & Ogiya (2018) for NFW haloes is typically a factor of 2 to 4 smaller than $r_{\rm conv}$ at fixed $N_{200}$. It is therefore unlikely to compromise the central mass profiles of simulated haloes. Furthermore, since $\epsilon_{\rm opt}\propto N_{200}^{-1/3}$ the ratio $\epsilon_{\rm opt}/m_{\rm DM}^{1/3}$ is fixed and $\epsilon_{\rm opt}$ can be conveniently expressed in units of the mean inter-particle distance: $\epsilon_{\rm opt}/l\approx 0.017$ (this is a factor $\approx 1.6$ larger than the maximum physical softening length for our fiducial runs). We suggest that softening lengths of order $\epsilon_{\rm opt}$ be employed in future large-scale simulations. This is a factor of 3 smaller than the radius for which we find better than 10 per cent convergence in circular velocity (eq. 20) and it is comparable to values adopted for most recent cosmological runs (Section 3.1). It is important to note, however, that this recommendation is based on our systematic study of convergence in the central mass profiles of DM haloes; other simulation statistics may prefer different $\epsilon$. 4. 4.

Convergence of halo mass functions: The mass functions of central CDM haloes are a robust prediction of cosmological simulations once a halo mass definition has been specified. For our fiducial softening lengths, Friends-of-friends ($M_{\rm FoF}$) and spherical overdensity ($M_{200}$) mass functions, converge to within $\approx 10$ per cent of those obtained from a higher-resolution run, above a mass scale corresponding to $\approx 32$ particles (Figure 13). Convergence is better at higher masses, reaching $\approx 5$ per cent for haloes resolved with at least 100 particles. These results are valid for a wide range of softening lengths, as demonstrated in Figure 13, where we compared the mass functions in a suite of 11 simulations that varied $\epsilon$ from $\approx 5.5\,{\rm pc}$ ($1/64^{th}$ of $\epsilon_{\rm fid}$) to $\approx 5.6\,{\rm kpc}$ (16 times $\epsilon_{\rm fid}$). For all but the largest softening length, all of these runs converge to within $\approx 5$ per cent for FoF mass, and to with $\approx 10$ per cent for $M_{200}$ provided haloes are resolved with at least 100 particles. Although small, the deviations depend systematically on $\epsilon$, particularly for haloes resolved with low numbers of particles (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100). For both mass definitions, the total abundance of haloes containing of order a few dozen to a few hundred particles increases systematically with increasing $\epsilon$, reaching a maximum at 4 to 8 times the fiducial softening (corresponding to roughly 4 to 8 per cent of the mean inter-particle separation) before rapidly declining (Figure 14). We emphasize, however, that such large softening lengths are undesirable as they compromise the innermost structure of DM haloes, limiting the simulation’s spatial resolution.

Overall, our results confirm and extend prior work on the accuracy and reliability of cosmological simulations of collisionless cold dark matter (e.g. P03; Navarro et al., 2004; Diemand et al., 2004; Lukić et al., 2007; Springel et al., 2008; Stadel et al., 2008; Power et al., 2016; van den Bosch & Ogiya, 2018; van den Bosch et al., 2018). The primary catalyst of this work was a follow-up study in which we repeat a number of these simulations but including either adiabatic hydro-dynamics, or a fully-calibrated set of sub-grid physics models for galaxy formation. In our opinion, it is therefore necessary to first highlight and expound the level of convergence in these pure DM runs, validating the numerical results in order to facilitate the interpretation of more complex hydrodynamical simulations.

We close by reiterating that our convergence study focused exclusively on “main” DM haloes, and that our results are unlikely to apply directly to sub-structure. As emphasized recently by van den Bosch & Ogiya (2018) and van den Bosch et al. (2018), the evolution of substructure in haloes is governed by several phenomena, both physical and numerical, that do not necessarily apply to isolated systems. We do, however, encourage future studies such as ours targeted explicitly at the structural evolution of substructure in hierarchical cosmologies.

Acknowledgments

We thank Lydia Heck and the cosma support team for their invaluable assistance, Alejandro Alejandro Benítez-Llambay and Matthieu Schaller for helpful conversations, and our anonymous referee for a thoughtful report. Our work has also benefited from various public python packages, including: scipy (Jones et al., 2001), numpy (van der Walt et al., 2011), matplotlib (Hunter, 2007) and ipython (Pérez & Granger, 2007). ADL acknowledges financial supported from the Australian Research Council through their Future Fellowship scheme (project number FT160100250). This work was supported by the Science and Technology Facilities Council [ST/P000541/1]. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk. This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.