Simulations of core formation for frequent dark matter self-interactions

TL;DR

This paper introduces a novel simulation method that models frequent dark matter self-interactions using adapted smoothed particle hydrodynamics, enabling more accurate studies of core formation in dark matter haloes.

Contribution

The authors develop a new $N$-body simulation framework that captures frequent dark matter self-interactions through derived energy transfer equations and thermal conductivity.

Findings

Simulated isothermal core formation in dark matter haloes.

Determined core growth rate as a function of self-scattering cross section.

Studied effects of different smoothing kernels.

Abstract

We present the first $N$-body simulations that adapt the equations of smoothed particle hydrodynamics to capture the effect of dark matter self-interactions which are too frequent to be resolved explicitly. The relevant energy transfer equations are derived, the appropriate thermal conductivity is determined and the effects of different smoothing kernels are studied. We apply our framework to simulate the formation of isothermal cores in isolated dark matter haloes and determine the core growth rate as a function of the self-scattering cross section. Our approach may be combined with explicit simulations of rare scatterings in order to simulate accurately the effects of arbitrary dark matter self-interactions in future cosmological simulations.

Equations

$$\frac{\Delta N}{N}=1-\exp \left(-\frac{\Sigma \sigma f}{m_{\text{DM}}}\right),$$

$$\frac{F_{\text{drag}}}{m_{\text{DM}}}\propto \frac{\sigma }{m_{\text{DM}}}\rho v_0^2.$$

$$\mathbf{j}=-\kappa \nabla T,$$

$$\rho \frac{\text{d}u}{\text{d}t}=-\nabla \mathbf{j},$$

$$\rho \frac{\text{d}u}{\text{d}t}=\nabla (\kappa \nabla T),$$

$$W(r,H)=\frac{8\pi }{H^3}\{\begin{array}{ll}1-6{\left(\frac{r}{H}\right)}^2+6{\left(\frac{r}{H}\right)}^3,& 0⩽\frac{r}{H}⩽\frac{1}{2},\\ 2{\left(1-\frac{r}{H}\right)}^3,& \frac{1}{2}<\frac{r}{H}⩽1,\\ 0,& \frac{r}{H}>1.\end{array}$$

$${\rho }_i=\sum _j{m}_{j}W(|{\mathbf{r}}_{ij}|,H_i).$$

$$\overline{A_i}=\frac{1}{{\rho }_i}\sum _j{A}_i{m}_{j}W(|{\mathbf{r}}_{ij}|,H_i).$$

$${\nu }_{x,i}^2=\overline{\Delta v_{x,i}^2}-{\overline{\Delta v_{x,i}}}^2,$$

$$\frac{\text{d}{u}_i}{\text{d}t}=\sum _j\frac{m_j}{{\rho }_i{\rho }_j}\frac{({\kappa }_j+{\kappa }_i)(T_j-T_i)}{|{\mathbf{x}}_{ij}{|}^2}{\mathbf{x}}_{ij}{\nabla }_i{W}_{ij},$$

$$\frac{\text{d}{u}_i}{\text{d}t}=2(\gamma -1)\sum _j\frac{m_j}{{\rho }_i{\rho }_j}\frac{{\kappa }_{ij}^{\prime }(u_j-u_i)}{|{\mathbf{x}}_{ij}{|}^2}{\mathbf{x}}_{ij}{\nabla }_i{W}_{ij},$$

$$\frac{\text{d}({\nu }_i^2)}{\text{d}t}=2(\gamma -1)\sum _j\frac{m_j}{{\rho }_i{\rho }_j}\frac{{\kappa }_{ij}^{\prime }({\nu }_j^2-{\nu }_i^2)}{|{\mathbf{x}}_{ij}{|}^2}{\mathbf{x}}_{ij}{\nabla }_i{W}_{ij}.$$

$$\frac{{\mathbf{v}}_i^{\prime }}{|{\mathbf{v}}_i^{\prime }|}=\frac{{\mathbf{v}}_i}{|{\mathbf{v}}_i|}.$$

$$\Delta {\mathbf{v}}_i=\left(\sqrt{\frac{v_i^{\prime 2}}{v_i^2}}-1\right){\mathbf{v}}_i.$$

$${\sigma }_{\text{V}}=\int {\sin }^2\theta \frac{\text{d}\sigma }{\text{d}\Omega }\text{d}\Omega ,$$

$${\lambda }_{\text{H}}=\frac{1}{n{\sigma }_{\text{V}}}=\frac{m_{\text{DM}}}{\rho {\sigma }_{\text{V}}},$$

$$\kappa \sim \frac{3k_B\nu }{2{\sigma }_{\text{V}}}$$

$$\frac{L}{4\pi r^2}=-\kappa \frac{\partial T}{\partial r}=-{\kappa }^{\prime }\frac{\partial {\nu }^2}{\partial r},$$

$$\frac{L}{4\pi r^2}=-\frac{3}{2}\frac{b}{\hat{a}}\rho \frac{{\lambda }_{\text{H}}^2}{t_r}\frac{\partial {\nu }^2}{\partial r},$$

$$\frac{L}{4\pi r^2}=-\frac{3b\nu }{2}{\left(\frac{{\sigma }_{\text{V}}}{m_{\text{DM}}}\right)}^{-1}\frac{\partial {\nu }^2}{\partial r},$$

$$\frac{L}{4\pi r^2}$$

$$=-\frac{3}{2}\hat{a}C\left(\frac{{\sigma }_{\text{V}}}{m_{\text{DM}}}\right)\frac{{\nu }^3\rho }{4\pi G}\frac{\partial {\nu }^2}{\partial r}.$$

$${\kappa }^{\prime }={\left(\frac{1}{{\kappa }_{\text{lmfp}}^{\prime }}+\frac{1}{{\kappa }_{\text{smfp}}^{\prime }}\right)}^{-1},$$

$${\kappa }_i^{\prime }=\frac{3}{2}\hat{a}b{\nu }_{r,i}\left(\frac{{\sigma }_{\text{V}}}{m_{\text{DM}}}\right){\left(\hat{a}{\left(\frac{{\sigma }_{\text{V}}}{m_{\text{DM}}}\right)}^2+\frac{b}{C}\frac{4\pi G}{{\rho }_i{\nu }_{r,i}^2}\right)}^{-1}.$$

$$M(r)=$$

$$\rho (r)=$$

$${\nu }_r^2(r)=$$

$$\displaystyle-\frac{r}{r+a}\left[25+52\frac{r}{a}+42\left(\frac{r}{a}\right)^{2}+12\left(\frac{r}{a}\right)^{3}\right]\Bigg{)}\;,$$

$$t_{\text{dyn}}=4\pi \sqrt{\frac{a^3}{GM}},$$

$$\hat{\sigma }$$

$$=0.05\left(\frac{{\sigma }_{\text{V}}/m_{\text{DM}}}{1{\text{cm}}^2{\text{g}}^{-1}}\right)\left(\frac{M}{10^{10}{\text{M}}_{\odot }}\right){\left(\frac{6.4\text{kpc}}{a}\right)}^2$$

$$t_0^{-1}$$

$${\chi }^2=\sum _j\frac{1}{{\overline{\sigma }}_j^2}{\left[{\overline{\text{d}({\nu }_r^2)/\text{d}t}}_j-{\left(\text{d}({\nu }_r^2)/\text{d}t\right)}_{\text{th},j}\right]}^2,$$

Full text

††pubyear: 2019††pagerange: Simulations of core formation for frequent dark matter self-interactions ††thanks: Preprint numbers: DESY-19-020, TTK-19-06–B

1 Introduction

The $\Lambda$CDM model based on the standard paradigm of collisionless cold dark matter (CDM) predicts the large-scale structure of the Universe remarkably well. But on small scales predictions from simulations of cosmological structure formation and observations appear to disagree with each other. In particular, pure CDM simulations predict cuspy density profiles with a steep slope ($\rho\propto r^{-1}$) in the centre of dark matter (DM) haloes (Dubinski & Carlberg, 1991; Navarro et al., 1996b, 1997), while observations tend to prefer more shallow profiles ($\rho\propto r^{0}$), corresponding to a central core (Flores & Primack, 1994; Moore, 1994; Moore et al., 1999a). This disagreement is known as the core-cusp problem. The observed number of subhaloes in the Milky Way is by far smaller than predicted by simulations (Moore et al., 1999b; Klypin et al., 1999). This discrepancy is known as the missing satellites problem. Another disagreement called the too-big-to-fail problem states that the most massive subhaloes in CDM simulations are too dense in the centre to host the observed satellites of the Milky Way (Boylan-Kolchin et al., 2011, 2012). Moreover, CDM is not able to explain the diversity of observed rotation curves of DM haloes with the same maximum circular velocity (Oman et al., 2015; Kamada et al., 2017). In addition to the other problems, this diversity problem constitutes the small-scale crisis of CDM (see the recent reviews by Bullock & Boylan-Kolchin (2017) and Tulin & Yu (2018)).

It has been suggested that these problems may be solved within the CDM framework by including baryonic processes such as star formation, supernova feedback and cooling into the $N$-body simulations (Navarro et al., 1996a; Governato et al., 2010), but whether this is sufficient in particular to solve the core-cusp problem is an ongoing debate (Read et al., 2018; Bose et al., 2018). As a result, modifications of the assumed properties of DM have been considered. Self-interacting dark matter (SIDM) was first suggested by Spergel & Steinhardt (2000) as another possibility to address the small-scale issues. The typical momentum-transfer cross section needed to solve these problems is of the order of magnitude $\sigma/m_{\text{DM}}\sim 1\ \text{cm}^{2}\text{g}^{-1}$.

At the same time bounds on the self-interaction cross section have been derived from a variety of different astrophysical systems, with the strongest ones arising from detailed studies of galaxy clusters. Bounds resulting from the non-observation of offsets between DM and galaxy centroids as strong as $\sigma/m_{\text{DM}}\lesssim 0.47~{}\text{cm}^{2}\text{g}^{-1}$ have been claimed (Markevitch et al., 2004; Randall et al., 2008; Kahlhoefer et al., 2014; Harvey et al., 2015). Subsequently, these bounds have been shown to be weaker than initially estimated, $\sigma/m_{\text{DM}}\lesssim 2~{}\text{cm}^{2}\text{g}^{-1}$  (Wittman et al., 2018) and are expected to be further diminished when taking into account the collisionality of galaxies (Kummer et al., 2018). Stronger bounds are obtained when considering sub-halo evaporation (Markevitch et al., 2004; Randall et al., 2008), $\sigma/m_{\text{DM}}\lesssim 1~{}\text{cm}^{2}\text{g}^{-1}$. Based on core-sizes of galaxy clusters bounds as strong as $\sigma/m_{\text{DM}}\lesssim 0.1~{}\text{cm}^{2}\text{g}^{-1}$ have been claimed, but these bounds strongly depend on assumptions of the baryonic physics (Kaplinghat et al., 2016; Elbert et al., 2018).

These results have all made the common assumption that a contact interaction between DM particles induces isotropic scattering. In this case the total cross section and the momentum-transfer cross section are very similar. Such models are relatively easy to implement in numerical studies, since self-interactions of this type have to be rare in order to satisfy observational constraints (Rocha et al., 2013; Peter et al., 2013). However, in order to address the small-scale problems and satisfy constraints from larger systems, it may be necessary to consider velocity-dependent self-interaction cross sections Buckley & Fox (2010); Loeb & Weiner (2011); Kaplinghat et al. (2016), such that the scattering rate increases with decreasing velocity. Such a velocity dependence is in fact well-motivated from particle physics and arises for example in models with light mediators Buckley & Fox (2010); Feng et al. (2010); Loeb & Weiner (2011); Tulin et al. (2013a, b).

While the velocity-dependence itself can be included in numerical simulations in a straight-forward way (Vogelsberger et al., 2012a; Zavala et al., 2013), it is no longer justified in these models to treat self-scattering as isotropic. In fact, for the case of a light mediator the differential scattering cross section is strongly peaked in the forward direction (Kahlhoefer et al., 2014). To implement such a model within the framework of rare self-interactions, it becomes necessary to introduce a cutoff angle and neglect the effect of small-angle scattering (Robertson et al., 2017b).

In this paper, we take a closer look at the effect of frequent small-angle scattering, which may give rise to a drag force and lead to the transport of energy (Yoshida et al., 2000; Moore et al., 2000; Gnedin & Ostriker, 2001; Markevitch et al., 2004; Ackerman et al., 2009; Feng et al., 2009). We develop a hybrid description for simulations of isolated haloes in the limit of frequent self-interactions. Using a hydrodynamical description we treat the DM particles as a collisionless fluid and implement the energy transport by solving a heat diffusion equation. We find that efficient heat conduction leads to a reduction of the central density of a DM halo, which can address the core-cusp problem. This approach is similar to the gravothermal fluid formalism that was developed to study globular clusters (Lynden-Bell & Eggleton, 1980) and has subsequently also been applied to self-interacting DM (Balberg et al., 2002; Ahn & Shapiro, 2005; Koda & Shapiro, 2011; Pollack et al., 2015; Essig et al., 2018; Choquette et al., 2018; Nishikawa et al., 2019).

We present an implementation of such a hybrid description for the energy transfer due to frequent self-interactions in the $N$-body code gadget-2 (Springel, 2005) using its formulation of smoothed particle hydrodynamics (SPH). As an application of this implementation, we study the evolution of cuspy DM haloes modeled with a Hernquist profile under the influence of frequent DM self-interactions. While we limit ourselves to systems with spherical symmetry in the present work, the framework can be generalised to non-symmetrical systems. Our implementation can be combined with existing implementations of rare self-interactions to perform simulations of SIDM that are valid for both large and small angle scatterings. This approach will enable robust predictions from simulations that will help to interpret observations more reliably.

This paper is organised as follows: In Sec. 2, we introduce frequent DM self-interactions. The implementation of energy transfer via frequent DM self-interactions is described in Sec. 3. The setup and initial conditions of our simulations are discussed in Sec. 4. In Sec. 5, we validate our simulations by comparing to analytic predictions. The results of our simulations are presented in Sec. 6 before we conclude in Sec. 7.

2 Frequent dark matter self-interactions

The strongest bound on the DM self-interaction cross section for velocity-independent scattering originates from the Bullet Cluster, more specifically from the requirement that the smaller DM halo (the subcluster) should not lose more that $20\texttt{--}30\%$ of its mass during core passage (Markevitch et al., 2004; Randall et al., 2008). Following Kahlhoefer et al. (2014), this statement can be expressed in terms of the number of scattering events leading to the expulsion of a DM particle from the subcluster (so-called expulsive collisions). In this context, a collision between a DM particle from the main cluster with a particle from the subcluster is expulsive if the velocity of both outgoing particles exceeds the escape velocity $v_{\text{esc}}$ of the subcluster.

For a subcluster moving through a larger cluster with relative velocity $v_{0}$ the fraction of expulsive collisions is found to be $f=(1-v_{\text{esc}}^{2}/v_{0}^{2})/(1+v_{\text{esc}}^{2}/v_{0}^{2})$ for the case of isotropic scattering. The fractional loss of DM particles can then be expressed as

[TABLE]

where $\Sigma$ is the integrated background density and $\sigma/m_{\text{DM}}$ the self-interaction cross section over DM mass.

In this approach, non-expulsive collisions are assumed not to affect the subcluster at all. This is a reasonable approximation if self-interactions are rare and $f$ is of order unity – a case that has been explored by a number of recent simulational studies (Robertson et al., 2017a; Kim et al., 2017). However, it is also possible to satisfy the bound from the Bullet Cluster for very large self-interaction cross sections provided that $f$ is sufficiently small. This is the case for example if the differential cross section is strongly peaked in the forward direction, as expected for interactions mediated by a light or massless exchange particle. In such a situation the majority of interactions do not lead to a significant transfer of energy and momentum, so that the total cross section (and hence the frequency of scattering processes) can be large without violating observational constraints.

For frequent self-interactions the majority of scattering processes do not lead to the immediate evaporation of DM particles, but they may affect the structure of the DM halo in different ways. To see this, consider a single DM particle, moving through a constant background density with velocity $\mathbf{v}_{0}$. After each scattering process with angle $\theta\ll 1$, the velocity in the direction parallel to $\mathbf{v}_{0}$ will be reduced by $\delta v_{\parallel}\approx v_{0}\sin^{2}\theta$, whereas in the direction perpendicular to $\mathbf{v}_{0}$ the velocity will increase by $\delta\mathbf{v}_{\perp}\approx v_{0}\sin\theta$. Over a large number of random scatters, $\delta\mathbf{v}_{\perp}$ will average to zero, but $\delta v_{\parallel}$ and $\delta v_{\perp}^{2}$ do not.

If we replace the single DM particle by an entire DM halo moving with velocity $v_{0}$, frequent self-interactions can therefore be thought of as causing a drag force:

[TABLE]

This drag force converts the directed motion of the DM halo into internal energy, i.e. random motion of the individual DM particles. Effectively, the DM halo is heated up, which can lead to the loss of particles via cumulative evaporation. This net effect of a large number of small-angle scatters is not captured when only considering the effect of rare scatterings with large momentum transfer.

In contrast to rare self-interactions, frequent self-interactions are notoriously difficult to implement in numerical simulations. The reason is that particles are typically expected to scatter multiple times during a single time-step, making an effective description of these interactions necessary. Indeed, since the mean free path $\lambda=(n\sigma)^{-1}$ is typically small ($\lambda\ll r$) for frequent self-interactions, we can approximately apply a fluid description. The idea is to model the energy transport with a heat diffusion equation for a pressureless fluid:

[TABLE]

where $\mathbf{j}$ denotes the heat flux and $\kappa$ the thermal conductivity. The temperature $T$ is assumed to be given by the one-dimensional velocity dispersion $\nu$ as $T=m_{\text{DM}}\nu^{2}/k_{B}$ (see e.g. Gnedin & Ostriker (2001)). Combined with energy conservation,

[TABLE]

the diffusion equation takes the form

[TABLE]

where $u$ is the internal energy per unit mass. As we will discuss below, the particle physics enters in the definition of the thermal conductivity $\kappa$, which depends in particular on the self-scattering cross section.

The idea to treat DM as a fluid has been used previously in the literature to study the structure (Gnedin & Ostriker, 2001; Yoshida et al., 2000; Moore et al., 2000) and stability (Balberg et al., 2002; Ahn & Shapiro, 2005; Koda & Shapiro, 2011; Pollack et al., 2015; Essig et al., 2018) of DM haloes for frequent self-interactions. The disadvantage of this approach compared to $N$-body simulations is that it cannot be directly extended to cosmological simulations and that one cannot easily include the effect of rare self-interactions with large momentum transfer, for which the effective description does not apply.

The present paper takes the first steps towards the inclusion of an effective description of frequent self-interactions in $N$-body simulations. We implement the resulting energy transfer in an existing code and validate our approach for a number of simple configurations. We discuss the scaling behaviour of the results and point out possible convergence issues. Our approach can then be combined with existing implementations of rare self-interactions, with the relative importance of the two effects depending on the underlying particle physics properties of the self-interaction cross section and the resolution of the simulation.

3 Implementation of energy transfer

In this section we present our implementation of frequent DM self-interactions in the cosmological $N$-body code gadget-2 (Springel, 2005). We adopt a hybrid description, in which DM is modeled as a collisionless fluid supplemented by an effective treatment of self-interactions. To calculate the energy transfer we use the SPH section of the code, which in fluid simulations is intended for computing the hydrodynamical forces. The central idea is to use a smoothing kernel $W(r,H)$ with kernel-support radius $H$ (i.e. the largest $r$ for which $W(r)>0$) to calculate local averages of relevant quantities in a way that achieves high resolution and low numerical noise. For concreteness, the cubic spline kernel is given by (Monaghan & Lattanzio, 1985)

[TABLE]

A comparison of different kernel choices is provided in App. A.

For a given smoothing kernel the local density of DM particles $\rho$ is calculated as

[TABLE]

We can then define mass- and kernel-averaged quantities as

[TABLE]

Using this approach we can calculate the one-dimensional velocity dispersions $\nu$:

[TABLE]

where $\Delta v_{x,i}=v_{x,i}-v_{x,j}$ is the relative velocity between the particle and its neighbouring particles in the kernel in the direction of interest. We calculate the velocity dispersion for the $x$-, $y$- and $z$-direction separately and obtain the total (3D) velocity dispersion $\nu_{3D}$ by summing the individual contributions. Assuming isotropy within the kernel, the radial velocity dispersion is simply given by $\nu_{r}^{2}=(1/3)\nu_{3D}^{2}$. The noise of $\nu_{r,i}^{2}$ is reduced by calculating the kernel average $\overline{\nu_{r,i}^{2}}$.

Frequent DM self-interactions cause a transfer of energy, which can be interpreted as a transport of the velocity dispersion $\nu$. We include this effect by solving the heat conduction equation given in Eq. (5). As shown by Jubelgas et al. (2004), its discretised version for the implementation in SPH codes reads

[TABLE]

where $u=\nu^{2}/(\gamma-1)=k_{B}T/((\gamma-1)m_{\text{DM}})$ is the thermal energy per unit mass (the one-dimensional case was discussed before by Brookshaw (1985)). With the replacement of $(\kappa_{i}+\kappa_{j})/2$ by $\kappa_{ij}=(2\kappa_{i}\kappa_{j})/(\kappa_{i}+\kappa_{j})$ to ensure a continuous heat flux (Cleary & Monaghan, 1999), we obtain

[TABLE]

where $\kappa^{\prime}_{ij}=(m_{\text{DM}}/k_{B})\kappa_{ij}$. The diffusion of the one-dimensional velocity dispersion is then given by

[TABLE]

At each time step this diffusion equation is solved for the radial velocity dispersion in order to obtain $\Delta\nu_{r,i}^{2}$, which corresponds to the kinetic energy transferred via frequent self-scatterings. The result, which can be positive or negative, then needs to be added to the squared velocity of particle $i$, ensuring energy conservation. Therefore, the velocity after energy transfer is given by $v_{i}^{\prime 2}=v_{i}^{2}+\Delta\nu_{r,i}^{2}$. Since thermal conduction is an effective description for energy transfer due to many scatterings for which the momentum transfers average out, the direction of the particle should not change, such that

[TABLE]

Hence we obtain

[TABLE]

To effect such a change in velocity, we calculate the conduction acceleration $\mathbf{a}_{i}=\Delta\mathbf{v}_{i}/\Delta t$ for a given time-step $\Delta t$ and add it to the gravitational acceleration. As usual, we also add a limiter to ensure that $v^{\prime 2}_{i}$ can never be smaller than zero.

The remaining challenge is to obtain a prescription for calculating the thermal conductivity $\kappa$ in terms of the underlying physics, in particular the differential scattering cross section $\mathrm{d}\sigma/\mathrm{d}\Omega$. Usually $\kappa$ is calculated from the mean free path, which depends on the total cross section $\sigma$. However, a more appropriate quantity in the present context is the viscosity cross section

[TABLE]

as discussed e.g. by Present (1958); McDaniel (1964); Massey et al. (1974); Schultz et al. (2008) and also in the context of DM by Tulin et al. (2013b); Cline et al. (2014). The viscosity cross section accounts for the fact that perpendicular scattering ($\theta=\pi/2$) is most efficient for heat transfer in a DM halo and for altering its density profile, while both forward and backward scattering have a reduced effect.111We note that for long-range interactions, the viscosity cross section is still logarithmically dependent on the small-angle cut-off, illustrating the fact that frequent interactions with small scattering angles dominate the energy transfer.

In analogy to the momentum transfer mean free path $\lambda_{\text{M}}=(n\sigma_{\text{T}})^{-1}$ defined in McDaniel (1964), we therefore define the heat transfer mean free path

[TABLE]

where $n$ is the number density. For a monatomic gas we have $\kappa\sim c_{v}\lambda_{\text{H}}n\nu$, where $c_{v}=3k_{B}/2$ is the specific heat per particle. We hence obtain for the conductivity

[TABLE]

This treatment is supported by a detailed calculation in Present (1958) that shows that the viscosity $\eta$ and thus the thermal conductivity ($\kappa=5/2\eta c_{v}$ for monatomic gases) scale as $\eta\sim\sigma_{\text{V}}^{-1}$.

To obtain a more precise expression for $\kappa$ we can employ the gravothermal fluid formalism used to study thermal conduction from DM self-interactions (Balberg et al., 2002; Ahn & Shapiro, 2005; Koda & Shapiro, 2011; Pollack et al., 2015; Essig et al., 2018). This formalism makes use of the flux equation

[TABLE]

where $L$ is the heat radiated through a sphere of radius $r$. Thus, we can obtain $\kappa^{\prime}$ from the relation between $L$ and the radial derivative of the velocity dispersion, replacing the total cross section by the viscosity cross section (Boddy et al., 2016).

We treat the short and long mean free path limits following Lynden-Bell & Wood (1968) and Balberg et al. (2002). Let us first consider the case that the mean free path $\lambda_{\text{H}}$ is much smaller than the gravitational scale height, i.e. the Jeans length $\lambda_{\text{J}}$ given by $\lambda_{\text{J}}^{2}=\nu^{2}/(4\pi G\rho)$. In this short mean free path regime the flux is given by

[TABLE]

where $t_{r}$ is the relaxation time. The constant $\hat{a}=\sqrt{16/\pi}\approx 2.257$ is the collision rate of DM particles (assuming a Maxwellian distribution) and $b=25\sqrt{\pi}/32\approx 1.385$ is derived from Chapman-Enskog theory (Lifshitz & Pitaevskii, 1981). The relaxation time can be written as $t_{r}=(\hat{a}\sigma_{\text{V}}\rho\nu/m_{\text{DM}})^{-1}$, leading to

[TABLE]

which fixes the constant of proportionality in Eq. (17).

If the heat transfer mean free path is larger than the Jeans length, particles can complete several orbits between two scatters in which significant energy is transferred. In this case the important length scale for conduction is the Jeans length $\lambda_{\text{J}}$ rather than the mean free path. In this long mean free path regime, the flux is given by

[TABLE]

The constant $C$ appearing in this equation must be determined by $N$-body simulations (Koda & Shapiro, 2011). A value of $C\approx 0.75$ is suggested for an NFW profile, while for a Hernquist profile the suggested value is $C\approx 0.9$ ($C\approx 1$ according to Balberg et al. (2002)). A very recent study finds $C=0.6$ for an NFW profile (Essig et al., 2018). As we will consider a Hernquist profile, we adopt $C=0.9$, emphasizing that uncertainties in $C$ enter linearly in the calculation of $\kappa^{\prime}$ in the long mean free path regime.

To interpolate between the short and long mean free path regimes, we take

[TABLE]

which yields

[TABLE]

We demonstrate in Appendix B that our code is able to reproduce the analytic solution for a one-dimensional test problem.

4 Simulating core formation in dark matter haloes

As discussed above, DM self-interactions cause energy transfer within DM haloes from regions with high temperature (i.e. high velocity dispersion) to colder regions. This effect heats up central cusps (which are colder than their surroundings) and transforms them into isothermal cores. For the remainder of this paper we will study this core formation, i.e. the transition from a cuspy profile to a cored profile. We will determine how the core size depends on both the age of the system and the strength of the self-interactions and determine appropriate rescaling rules.

To initialise our simulations, we model individual DM haloes with cuspy Hernquist profiles (Hernquist, 1990), for which the enclosed mass, density profile and one-dimensional velocity dispersion are given by

[TABLE]

where $M$ is the total mass and $a$ is the scale radius. The $N$-body solver furthermore requires the introduction of a softening length, and we choose the Plummer equivalent softening length according to Power et al. (2003). We use the tool makegal to generate the initial conditions for spherical haloes in dynamical equilibrium (Springel & White, 1999). makegal first samples the mass profile using the analytic expression Eq. (25). The velocities of the DM particles are drawn from random samples to sample the distribution function of the Hernquist profile.

We focus on a DM halo with $M=10^{10}~{}\text{M}_{\odot}$, corresponding to the size of a typical dwarf galaxy. The scale radius $a=6.4~{}\text{kpc}$ is chosen to satisfy the concentration-mass relation from Dutton & Macciò (2014). We define the dynamical time as the orbital time at the scale radius as in Kochanek & White (2000); Robertson et al. (2017b):

[TABLE]

implying $t_{\text{dyn}}=0.98~{}\text{Gyr}$ for the halo that we consider. We run simulations with $N=10^{6}$ and $N=10^{7}$ particles. For the smaller number of particles, numerical cores, which form due to gravitational force softening (Robertson et al., 2017a), start to grow for $t>5~{}t_{\text{dyn}}$. For the case of $N=10^{7}$ stability of the initial conditions is guaranteed for the entire timescale we simulate. Adding the conduction module to the SPH version of gadget-2 increases the runtime of the code by roughly 30 per cent.

Following Koda & Shapiro (2011); Vogelsberger et al. (2012b) we also introduce a dimensionless cross section

[TABLE]

and a dimensionless time

[TABLE]

where $\rho_{0}=M/(2\pi a^{3})$ and $v_{0}=a\sqrt{4\pi G\rho_{0}}$.

We find that for cross-sections $\hat{\sigma}\lesssim 3$ the entire halo (down to $r\sim 0.01a$, which are the smallest radii that can be resolved in the simulation) resides in the long mean free path regime for the timescales we simulate. This implies that the thermal conductivity $\kappa$ is directly proportional to the self-interaction cross section and hence we expect larger self-interactions to lead to larger core sizes and more rapid core formation. Finally, we note that in the long mean free path limit $\kappa$ depends on the combination $C\sigma_{\text{V}}$, where $C$ is an empirical parameter with a considerable uncertainty (see Eq. (22) and surrounding text). Results for different values of $C$ can be obtained from the results presented below by an appropriate rescaling of $\sigma_{\text{V}}$.

5 Validation

To validate the implementation presented above we compare the heat fluxes calculated by the code immediately after initialisation (i.e. before applying the first time-step) to the analytic expectation for a Hernquist profile. Specifically, we compare the respective rates $\text{d}(\nu_{r}^{2})/\text{d}t$ obtained from the solution of the heat equation, see Eq. (12). For the numerical result we calculate the median of all particles in a radial shell to reduce noise. The comparison with the analytic result is shown in Fig. 1 for $\sigma_{\text{V}}/m_{\text{DM}}=1~{}\text{cm}^{2}\text{g}^{-1}$. In the left panel, we fix the total number of particles in the simulation to $N=10^{6}$ and consider different numbers of neighbours $N_{\text{nbg}}$ in the smoothing kernel, while in the right panel, both $N$ and $N_{\text{nbg}}$ are varied.222At first sight, the fact that $\mathrm{d}(\nu^{2})/\mathrm{d}t$ is negative in the central region is counter-intuitive, given that we expect an increase in temperature in this region. The reason for this behaviour is that a self-gravitating system has negative heat capacity. In other words, reducing the kinetic energy of particles in the centre leads to an increase in total energy due to the virial theorem. This forces the particles onto larger orbits and drives the transformation of a cusp into a core (Pollack et al., 2015).

We observe that for values of $N_{\text{nbg}}$ as typically chosen for SPH problems ($N_{\text{nbg}}\sim 10^{2}$) we substantially overestimate the conduction effect. The reason is that there are no additional hydrodynamical forces that reduce the noise in the velocity distribution. We thus need to increase the number of neighbours for the cubic spline kernel and thereby the kernel size to achieve a sufficient reduction in noise.333Note that we can allow a large number of neighbours even for the cubic spline kernel because we do not compute hydrodynamical forces and therefore are not limited by pairing instabilities (Dehnen & Aly, 2012). At the same time making $N_{\text{nbg}}$ too large partially averages out real physical effects. This reasoning implies the existence of an optimum range for $N_{\text{nbg}}$.

We estimate this optimum value using the least-squares method, i.e. by minimising

[TABLE]

where $\overline{\text{d}(\nu_{r}^{2})/\text{d}t}_{j}$ is the median of the numerical results in the $j$th radial shell and $(\text{d}(\nu_{r}^{2})/\text{d}t)_{\text{th},j}$ is the theoretical prediction for a Hernquist profile. The uncertainty $\overline{\sigma}_{j}$ of the median is difficult to determine, since for neighbouring particles the values of $\text{d}(\nu_{r}^{2})/\text{d}t$ are highly correlated. If $\sigma_{j}$ denotes the standard deviation of $\text{d}(\nu_{r}^{2})/\text{d}t$ within each radial shell, we expect $\overline{\sigma}_{j}=\epsilon\,\sigma_{j}/\sqrt{N_{\text{sh},j}}$, where $N_{\text{sh}}$ is the number of particles in the radial shell and $\epsilon>1$ is some correction factor that accounts for correlations. Since we will only be interested in determining the best-fit value of $N_{\text{nbg}}$ and not in a goodness-of-fit estimate, we can drop this factor and simply determine the minimum of

[TABLE]

For the fitting procedure we focus on a limited range around $r/a\sim 1$, which ensures on the one hand that there is a sufficiently large number of particles in each shell and on the other hand that the size of the smoothing kernel does not bias our results. For $N=10^{6}$ we achieve the best agreement with the analytic expectation in this range with $N_{\text{nbg}}=655$. Increasing the total number of particles $N$ in the simulation induces more noise in the distribution of $\text{d}(\nu_{r}^{2})/\text{d}t$. To compensate for this effect we must increase $N_{\text{nbg}}$ (see right panel of Fig. 1). We find $N_{\text{nbg}}$ to scale approximately proportional to $N^{0.4}$ and obtain an optimal vaule of $N_{\text{nbg}}=1700$ for $N=10^{7}$. We emphasize that numerical convergence in SPH is only possible in the simultaneous limit $N\rightarrow\infty$ and $N_{\text{nbg}}\rightarrow\infty$ (Zhu et al., 2015).

Moving to more extended kernels (such as Wendland C6 as proposed by Dehnen & Aly (2012) or Wendland C8) does not improve the agreement compared to the cubic spline kernel, but leads to higher computational cost. Details of the kernel comparison are presented in App. A.

For very small radii, $r\ll a$ some discrepancies between analytical and numerical result remain even with the optimum choice for $N_{\text{nbg}}$. The reason is that the requirement of a large number of neighbours implies a large kernel-support radius $H$. Once $H$ becomes comparable to the radial distance of the particle to the centre of the halo, our implementation no longer produces reliable results, because the particles sampled by the kernel no longer provide unbiased estimates of the underlying distributions. In other words, it may no longer be possible to find a sufficient number of particles with smaller radius to reliably estimate the gradient in the velocity dispersion.

Fortunately, the limited resolution in the central region does not impede our ability to study core formation. The reason is that the core formation rate depends on the total amount of heat transported towards the centre. This heat transport can be accurately estimated from our implementation, which agrees well with the analytic expectations for $r\gtrsim a/2$. However, the numerical artefacts in the central region may lead to unphysical energy transport within the core, which may affect core stability. We therefore cannot cover the full evolution of the halo in the gravothermal fluid model and do not expect to accurately capture the recollapse of cores at very late times (Lynden-Bell & Wood, 1968; Balberg et al., 2002; Koda & Shapiro, 2011; Pollack et al., 2015; Kochanek & White, 2000; Robertson et al., 2017b).

6 Results

Based on our findings in the previous section, we use the cubic spline kernel with $N_{\text{nbg}}=655$ for $N=10^{6}$ and $N_{\text{nbg}}=1700$ for $N=10^{7}$ to simulate the dwarf-size DM halo defined in Sec. 4. Fig. 2 shows the density profile and velocity dispersion of the halo after time $t_{0}$, defined as the relaxation time at the centre. In agreement with expectation we find heat conduction from frequent DM self-interactions to reduce the central density and to start core formation. At the same time the velocity dispersion flattens in the centre of the halo, leading to an isothermal profile. The fact that the velocity dispersion is not exactly flat in the centre stems from the numerical limitations discussed in Sec. 5. The results for $N=10^{6}$ and $N=10^{7}$ agree very well, indicating that our simulations converge for $\sigma_{\text{V}}/m_{\text{DM}}=1~{}\text{cm}^{2}\text{g}^{-1}$. On purpose, our fiducial choice of particles is not the maximal number that we can afford but $N=10^{6}$ and $N_{\text{nbg}}=655$ because our technique is intended to be used in larger simulation boxes where only moderate particle numbers per halo can be afforded. A comparison to simulations with $N=10^{7}$ particles supports the presented results with a deviation of the final core size smaller than 2 per cent for $\hat{\sigma}=1$.

Fig. 3 illustrates the time evolution of the DM halo in dimensionless units for $\hat{\sigma}=1$ after $t=2.6~{}t_{0}$, $t=12.8~{}t_{0}$ and $t=25.5~{}t_{0}$. We extract the core size of the DM halo by fitting the density profile to a cored Hernquist profile

[TABLE]

We keep the parameter $\beta=4$ fixed and allow the other parameters to vary during the fitting processes (Robertson et al., 2017b). The core growth for different dimensionless cross sections $\hat{\sigma}$ is shown in the left panel of Fig. 4. The growth of the core is very rapid in the beginning and decreases with time. For comparison, Pollack et al. (2015) found a maximal core size of a similar magnitude for an NFW profile as we find at $t\approx 25~{}t_{0}$ with the same cross section444Note the difference in the definition of the dimensionless cross section in Pollack et al. (2015) $\hat{\sigma}^{\prime}=\hat{\sigma}/(4\pi)$ such that $\hat{\sigma}^{\prime}=0.08$ is equivalent to $\hat{\sigma}=1$..

In the right panel of Fig. 4 we show the core size as a function of the product of dimensionless time and cross section. The fact that all curves lie on top of each other implies that the growth rate depends linearly on the cross section. This scaling originates from the fact that the entire halo is in the long mean free path limit (see discussion in the end of Secs. 3 and 4).

We find that the core size at $t=4~{}t_{\text{dyn}}$ is $r_{\text{core}}/a\approx 0.45$ for $\hat{\sigma}=1$. This is in agreement with the results of Kochanek & White (2000) as they find a maximal core size of $r_{\text{core}}/a\approx 0.4$ for the same cross section. This core size already reached at $t\approx t_{\text{dyn}}$, in our case the evolution is much slower. Note that the maximal core size found by Robertson et al. (2017b) is much smaller ($r_{\text{core}}/a\approx 0.15$) compared to our results and the results of Kochanek & White (2000).

Since we are not able to resolve the very centre of the halo, we are not able to simulate the full evolution of the gravothermal fluid model. Especially for large cross sections a recollapse of the core is expected for the timescale we simulate. For example, for $\hat{\sigma}=3$ we have $t_{0}\approx 0.05~{}\text{Gyr}$, which means we cover $t=4\,t_{\text{dyn}}\approx 80\,t_{0}$ in our simulation, while Vogelsberger et al. (2012b) found that the recollapse starts after $t\approx 25\,t_{0}$. The prediction of the core size for the large cross section $\hat{\sigma}=3$ is therefore not reliable, since we expect the core to recollapse. Note that $\hat{\sigma}=3$ corresponds to a very large cross section of $\sigma_{\text{V}}/m_{\text{DM}}\approx 59~{}\text{cm}^{2}\text{g}^{-1}$ which is not compatible with current bounds. For cross sections which are compatible with the bounds and for timescales smaller than the age of the Universe recollapse is not relevant for isolated DM haloes.

7 Discussion

We have presented the first $N$-body simulations using SPH to capture the effect of DM self-interactions for frequent scattering. Using an implementation of thermal conduction in the SPH section of gadget-2 we carried out simulations of an (isolated) dwarf-sized DM halo. To map the particle physics parameters to the thermodynamical quantities (i.e. the thermal conductivity) we use the well-established idea of the gravothermal fluid formalism. The resulting energy transfer from the outer part of the halo towards the centre has a significant impact on the evolution of the system. The central density of an initially cuspy halo is reduced and an isothermal core is created. We therefore recover the well-known mechanism for solving the core-cusp problem and the too-big-to-fail problem with SIDM.

Due to the negative heat capacity of self-gravitating systems, any self-interacting DM halo will experience a gravothermal catastrophe (Lynden-Bell & Wood, 1968) at very late times or for very large cross sections (Balberg et al., 2002; Kochanek & White, 2000). The fact that the resulting DM haloes would be in obvious conflict with observations can be used to obtain relevant constraints on the self-interaction cross section, in particular when including effects from dissipation (Essig et al., 2018). While it would be very interesting to investigate these effects in our simulations, numerical limitations in the central region of the DM halo imply that it is currently not possible to simulate its recollapse. It remains to be seen whether this problem can be overcome with higher-resolution simulations or different smoothing kernels.

Our formalism provides a proof of principle for the inclusion of frequent DM self-interactions in $N$-body simulations and should be seen as a first step towards a comprehensive treatment of all effects of SIDM. For example, our simulations do not presently include the effective drag force created by frequent self-interactions, i.e. we simulate the SIDM effects on the halo via heat conduction only. The reason for neglecting this drag force is that it is expected to be small in the system with spherical symmetry, as the velocity distribution of an isolated halo and therefore in each kernel is isotropic. In order to study anisotropic systems such as merging galaxy clusters, an implementation of an effective drag force will be needed.

Furthermore, the set-up of our simulations does not fully capture the relevant physical process during structure formation. Considering initially cuspy profiles implicitly assumes that self-interactions are not relevant for structure formation, which may not be a good approximation for large self-interaction cross sections. In order to extract reliable predictions and bounds from our formalism it will be essential to perform full cosmological simulations. The fact that such simulations involve many different length scales complicates the optimum choice of the smoothing kernel and may necessitate adaptive kernel sizes.

Finally, we have focussed on simulations of DM only. Understanding the interplay of baryonic physics and SIDM is an important next step, in particular in the context of the core-cusp problem. For simulations of rare DM self-interactions there has been a lot of progress on this subject recently. Testable predictions of SIDM have been developed and a preference for small self-interaction cross-sections ($\sigma/m_{\text{DM}}\approx 0.1~{}\text{cm}^{2}\text{g}^{-1}$) from core sizes in galaxy clusters has been found (Elbert et al., 2018; Robertson et al., 2018a; Sameie et al., 2018; Despali et al., 2018) although baryonic feedback introduces a sizeable scatter (Robertson et al., 2018b).

It is important to emphasize that our approach for simulating frequent self-interactions should not be seen as an alternative to existing state-of-the-art simulations of SIDM. Indeed, our ultimate goal is to combine the implementation presented here with the well-established formalism for rare scattering in order to simulate self-scattering cross-sections with arbitrary angular dependence. The idea would be to explicitly simulate large-angle scatters while using our approach to capture the effect of small-angle scatters. This hybrid solution will then make it possible to perform fully consistent $N$-body simulations of observationally favoured SIDM models with scattering cross sections depending on velocity and scattering angle. Doing so will enable us to draw robust conclusions on the nature of DM self-interactions.

Acknowledgements

We thank Manoj Kaplinghat and Denis Wittor for valuable discussions. This work is funded by the Deutsche Forschungsgemeinschaft (DFG) through the SFB grant 676, the Emmy Noether Grant No. KA 4662/1-1, the ERC Starting Grant ‘NewAve’ (638528) and under Germany’s Excellence Strategy – EXC 2121 ,,Quantum Universe“ – 390833306.

Appendix A Kernel comparison

We have experimented with different choices for the SPH kernel in our simulation. In particular, we have investigated the effect of more extended smoothing kernels, such as the Wendland C6 kernel (WC6) (Dehnen & Aly, 2012) defined by

[TABLE]

where $(\cdot)_{+}=\operatorname{max}\{0,\cdot\textperiodcentered\}$. More extended kernels obviously have the property of a larger effective kernel size (kernel-support radius $H$ i.e. the largest value of the radius, $r$, for which $W(r)>0$). This increases the number of neighbour particles, while keeping the same smoothing scale $h$ (Price, 2012). The smoothing scale is essentially defined by the kernel standard deviation $h=2\sigma$. A typical choice for the WC6 kernel is $N_{\text{nbg}}=295$. Finally, we have also experimented with a WC8 kernel (Wendland, 2009) defined as

[TABLE]

Keeping the smoothing scale, $h$, of the WC6 kernel with $N_{\text{nbg}}=295$ fixed, the WC8 kernel allows for $N_{\text{nbg}}=382$. We find that this choice leads to no significant improvement compared to the cubic spline kernel for the agreement of $\text{d}(\nu_{r}^{2})/\text{d}t$ with the theoretical expectation $(d(\nu_{r}^{2})/dt)_{\text{th}}$ as shown in Fig. 5. However, the computational cost increases substantially for the same smoothing length, $h$. Consequently, we need to increase the smoothing scale at the cost of losing resolution. This is achieved by increasing the number of neighbors for a fixed kernel. A useful quantity for the efficiency of the kernel is the number of neighbours inside the smoothing scale $N_{h}=(h/H)^{3}N_{\text{nbg}}$. $N_{h}$ is larger for the cubic spline kernel compared to the more extended kernels for a fixed number of particles. To keep the computational cost low we choose the cubic spline kernel for the simulations presented in Sec. 6. To summarise, we learned that the smoothing scale of the kernel needs to be increased by raising the number of neigbours $N_{\text{nbg}}$ for a given kernel in order to reduce the noise effectively.

Appendix B One-dimensional test problem

For the purpose of testing our implementation of heat conduction, we set up a one-dimensional test without gravity and compared the results to the analytical solution. To this end, we initialised particles on a grid with random velocities but different velocity dispersions on either side of the centre of the computational domain. We calculate the velocity dispersions from the initial velocities. The velocities are updated owing to the conduction and, hence, the velocity dispersions change accordingly. However, in this particular test we did not allow the particles to move but instead kept them fixed on the grid. This way, we can study only the transport of heat and compare the results to the analytic solution. Initially, we set the conductivity $\kappa^{\prime}$ to a certain value and keep it constant throughout the whole test $\kappa^{\prime}=\alpha=10$ since $\alpha=\kappa^{\prime}/\rho$ and we chose $\rho=1$ in internal units. Converted to physical units, this implies $\rho=6.77\cdot 10^{-22}~{}\text{g}/\text{cm}^{3}$, $\kappa^{\prime}=2.04\cdot 10^{6}~{}\text{g}/(\text{cm}\cdot\text{s})$ and $\alpha=3.02\cdot 10^{27}~{}\text{cm}^{2}/\text{s}$. For this case, the analytic solution is known (Jubelgas et al., 2004)

[TABLE]

where $x_{0}$ is the position of the initial difference of size $\Delta\nu$ of the velocity dispersion, and $\nu_{0}$ the mean velocity dispersion. We use the WC6 kernel with $N_{\text{nbg}}=295$ and find that the resulting thermal conduction of the simulation agrees well with the analytic solution around the centre for constant $\kappa^{\prime}$. The results are shown in Fig. 6.