Gravothermal collapse of isolated self-interacting dark matter haloes: N-body simulation versus the fluid model

Self-Interacting Dark Matter (SIDM) is a collisional form of cold dark matter (CDM), originally proposed to solve problems that arose when the collisionless CDM theory of structure formation was compared with observations of galaxies on small scales. The quantitative impact of the proposed elastic collisions on structure formation has been estimated previously by Monte Carlo N-body simulations and by a conducting fluid model, with apparently diverging results. To improve this situation, we make direct comparisons between new Monte Carlo N-body simulations and solutions of the conducting fluid model, for isolated SIDM haloes of fixed mass. This allows us to separate cleanly the effects of gravothermal relaxation from those of continuous mass accretion in an expanding background universe. When these two methods were previously applied to halo formation with cosmological boundary conditions, they disagreed by an order of magnitude about the size of the scattering cross section required to solve the so-called 'cusp-core problem.' We show here, however, that the methods agree with each other within 20 per cent for isolated haloes. This suggests that the two methods are consistent, and that their disagreement for cosmological haloes is not caused by a breakdown of their validity. The isolated haloes studied here undergo gravothermal collapse. We compare the solutions calculated by these two methods for gravothermal collapse starting from several initial conditions. This allows us to calibrate the heat conduction which accounts for the effect of elastic hard-sphere scattering in the fluid model. The amount of tuning of the thermal conductivity parameters required to bring the two methods into close agreement for isolated haloes, however, is too small to explain the discrepancy found previously in the cosmological context.


INTRODUCTION
In the currently standard cosmological model (ΛCDM), a flat universe with a cosmological constant contains collisionless cold dark matter as its dominant matter component, perturbed by primordial Gaussian-random-noise density fluctuations. This model has been highly successful at explaining observations of the background universe and large-scale structure. On small scales, however, the distribution of dark matter in coordinate and velocity space is not E-mail: junkoda@physics.utexas.edu. † E-mail: shapiro@astro.as.utexas.edu fully understood. N -body simulations show that the density profiles of the virialized regions ('haloes') that form in this collisionless dark matter are cuspy, such as the Navarro, Frenk, & White (1997, NFW) profile, in which ρ → r −1 toward the centre. Recent high-resolution simulations show that the inner profile is not exactly a power law (Hayashi et al. 2004;Diemand et al. 2005;Merritt et al. 2006;Gao et al. 2008;Stadel et al. 2009;Navarro et al. 2010), but it diverges nevertheless. On the other hand, a cored profile (a profile which flattens in the centre) such as a pseudoisothermal profile, is favored by observations of dwarf and low surface brightness (LSB) galaxies (de Blok et al. 2003;Gentile et al. 2004;Kuzio de Naray et al. 2006;Zackrisson et al. 2006;Gentile et al. 2007;Kuzio de Naray et al. 2008;de Blok et al. 2008;Oh et al. 2010b;de Blok 2010). Dwarf spiral and LSB galaxies are dark-matter dominated. As such, it was originally thought, their mass distribution should reflect the dark matter dynamics alone, and be relatively less affected by the complexity of the dissipative baryonic component. This, it was thought, makes these systems ideal for studying the undisturbed, intrinsic, dark matter distribution on small scales.
This apparent cusp-core conflict was one of the smallscale structure problems of the CDM model which prompted suggestions that the dark matter might be something else, with microscopic properties that would alter the structure on small scales without spoiling the success of CDM on largescales. Spergel & Steinhardt (2000) proposed self-interacting dark matter (SIDM) as a possible solution to this cusp-core problem, adding hypothetical elastic-scattering collisions to the otherwise collisionless particles of the standard CDM cosmology. Heat transfer within the virialized haloes, due to these non-gravitational collisions, was then suggested to make the halo cores expand. This latter idea was confirmed by several numerical and analytical studies. Burkert (2000) introduced a Monte Carlo scattering algorithm between dark matter particles to take the self interaction into account in a numerical N -body simulation, and this method was refined by Kochanek & White (2000). These N -body results suggested, however, that SIDM haloes would undergo gravothermal collapse, making them unsuitable to explain the observed haloes. Balberg, Shapiro, & Inagaki (2002, hereafter BSI) applied a conducting fluid model, originally invented to describe gravitational scattering in star clusters, to isolated 1 SIDM haloes. BSI derived a self-similar gravothermal collapse solution at large Knudsen number (Kn 1, where Kn is the ratio of SIDM scattering mean free path to the system size). The isolated halo collapsed within a finite time. They also showed that the collapse is delayed compared to their self-similar solution when the Knudsen number is comparable to or smaller than one, because the length scale of energy exchange is restricted by the mean free path. Their conclusion that the collapse time would naturally exceed a Hubble time was more optimistic about the SIDM hypothesis.
A realistic halo is not isolated, however, since it forms in a cosmologically expanding background universe, with infall and a finite pressure at the virial radius (Shapiro et al. 1999). Cosmological simulations showed that a cross section per unit mass σ = 0.5 − 5 cm 2 g −1 makes the profile cored, and the cored profile is stable (Yoshida et al. 2000;Davé et al. 2001). Colín et al. (2002) emphasized that the profile depends on the accretion history, especially when the last major merger occurred.
With the importance of cosmological infall in mind, Ahn & Shapiro (2005) derived a cosmological similarity solution for this problem, with proper account taken of such cosmo-logical boundary conditions. This solution shows that the gravothermal collapse, which occurs for the isolated halo, is prevented by infall in the cosmological case, and the core has a constant size in units of the virial radius, for a given SIDM cross section. When there is no self-interaction, this fluid approximation gives a density profile similar to the cuspy profile found in N -body simulations. This shows that the fluid approximation also describes the virialization of collisionless dark matter appropriately, providing in effect an analytical derivation of the NFW profile in the collisionless limit. This is because the collisionless Boltzmann equation reduces to fluid conservation equations for an ideal gas with ratio of specific heats equal to 5/3 if the velocity distribution of the particles in the frame of bulk motion is isotropic and skewless (see Section 2.2). In the presence of SIDM scattering collisions, too, those analytical solutions are in qualitative agreement with the corresponding Monte Carlo N -body simulations; SIDM haloes have cores and the cores collapse within a finite time when they are isolated, but they do not collapse within a Hubble time in a cosmological environment. However, the values of the cross section necessary to explain the observed dark matter density profiles are not in agreement. Ahn & Shapiro (2005) find that σ ∼ 200 cm 2 g −1 fits the dwarf and LSB galaxy rotation curves best, while N -body simulations suggest the range of σ = 0.5 − 5 cm 2 g −1 gives the observed central density.
In order to test the SIDM hypothesis by comparison with astronomical observations, it is necessary to improve our understanding of these theoretical predictions and reconcile them. That will be the focus of this paper, as described below. In the meantime, related progress continues to be made on other fronts, in comparing both SIDM and collisionless CDM with observation. For example, non-circular motions may affect the density profile estimate (Hayashi & Navarro 2006), but they are usually not strong enough to make observations consistent with the theoretically-predicted cuspy profile Trachternach et al. 2008;Kuzio de Naray et al. 2009). The cusp-core conflict may be mitigated by baryonic processes of gas outflow, induced by supernovae feedback (Governato et al. 2010;Oh et al. 2010a). Central dark-matter densities could also be reduced, in addition, by bars (Weinberg & Katz 2002), or by gas clumps sinking via dynamical friction (El-Zant et al. 2001), while the system is baryon rich, but these two scenarios may not be strong enough to make large enough cores in realistic situations (Sellwood 2008;Jardel & Sellwood 2009). While observations of dwarf and LSB galaxies generally prefer cored density profiles, some, however, are also consistent with an NFW profile (Hayashi et al. 2004;Simon et al. 2005;Hayashi & Navarro 2006;Valenzuela et al. 2007). The wide diversity of dwarf/LSB cores has also led to a suggestion that SIDM alone cannot be the full solution to the cusp-core problem (Sánchez-Salcedo 2005;Kuzio de Naray et al. 2010).
There are also some constraints on the value of the SIDM cross section from comparison of theoretical predictions with galaxy clusters, for which dark matter velocity is much higher than for galaxies. N -body results find that the core size of relaxed SIDM clusters of galaxies becomes too large for cross section σ 1 cm 2 g −1 (Yoshida et al. 2000;Arabadjis et al. 2002;Lewis et al. 2003). The analytical cosmological self-similar solutions of Ahn & Shapiro (2005), however, point out that the core sizes of clusters are also small enough, not only for small cross section in the long-mean-free-path regime, but also for large cross section in the short-mean-free-math regime, i.e., σ 200 cm 2 g −1 , because the short mean free path then limits the amount of heat conduction. However, such a large cross section would also enhance the fluid-like behaviour of SIDM, which may then conflict with observations of merging clusters. For example, observations of cluster 1E 0657-56, known as the 'bullet cluster,' from which total matter density has been mapped by weak and strong gravitational lensing measurements while the density of the intergalactic baryon-electron fluid was mapped by X-ray measurements, show that the dark matter spatially segregated from the baryon-electron plasma, as it would be if the dark matter is not highly collisional (Clowe et al. 2006). The dark matter and galaxies of the subcluster have passed through the main cluster without distortion while the baryon gas shows a bow shock due to its collisional nature. This observation excludes the possibility that SIDM is too highly collisional, or fluid-like. Analytical estimates and Monte Carlo N -body simulations of the bullet cluster by Markevitch et al. (2004) and Randall et al. (2008) constrain the velocity independent cross section to be σ < 0.7 cm 2 g −1 , using the fact that the massto-light ratio of the subcluster is normal. Larger cross section would make the mass-to-light ratio smaller because the SIDM collisions scatter the dark matter out of the subhalo.
[On the other hand, see Mahdavi et al. (2007) for a possibility that SIDM with cross section ∼ 4 cm 2 g −1 explains a cluster Abell 520 which has substructures with anomalous mass-to-light ratio]. If the cross section is velocity dependent, it puts a constraint only at very large relative velocity. The relative velocity of the merging haloes is estimated to be between 2, 500 and 4, 000km s −1 at their observed separation, and even higher at centre passage (Milosavljević et al. 2007;Springel & Farrar 2007;Mastropietro & Burkert 2008). If the cross section is too large, the SIDM halo could also be too spherical compared to the observed elliptical matter distribution of galaxy clusters (Miralda-Escudé 2002). In addition, large cross sections result in too much heating and evaporation of substructure haloes in clusters (Gnedin & Ostriker 2001). These constraints almost exclude the possibility that a velocity-independent cross section solves the cusp-core problem for dwarf/LSB galaxies, but, velocity-dependent cross section can still be effective on dwarf scales, and simultaneously negligible on galaxy cluster scales: e.g., inversely proportional to relative speed (Firmani et al. 2000;, or short range Yukawa interaction, whose scattering cross section has v −4 velocity dependence for large relative velocities (Koda 2009;Loeb & Weiner 2010).
There was an additional motivation for the SIDM hypothesis when it was first put forth, involving the overabundance of subhaloes in CDM N -body simulation compared to observations of the Local Group. The number of substructures in the collisionless CDM model has previously been thought to be about an order of magnitude larger than the observed number of dwarf galaxies in the Local Group. Although the number of substructures is reduced by SIDM stripping, D'Onghia & Burkert (2003) claim that SIDM does not solve this problem; the cross-section that makes the profile cored (0.6 cm 2 g −1 ) is apparently not efficient enough to reduce the number of satellites down to the observed level. Recent findings of nearby faint galaxies improve the agreement between observations and collisionless CDM simulations, but still require some feedback or reionization effect that suppresses star formation in low-mass haloes (e.g. Koposov et al. 2009;Wadepuhl & Springel 2010, and references therein).
Apart from the original motivation for proposing SIDM as a solution of small-scale structure problem of the collisionless CDM model, interest in the SIDM hypothesis remains strong for other reasons, as well. The fundamental nature of dark matter is still unknown. Anomalous cosmic ray detections, reported by PAMELA, FERMI and ATIC, have recently stimulated new particle physics models for dark matter that result in the scattering interaction of SIDM as a secondary effect (Arkani-Hamed et al. 2009;Buckley & Fox 2010;Feng et al. 2010). Astronomical constraints on the SIDM cross section and its velocity dependence continue to be of importance, therefore, in order to constrain such models.
Further progress along these lines requires that we reduce the uncertainties in the theoretical productions for SIDM. As described above, the quantitative estimates by Monte Carlo N -body simulation and the conducting fluid model, of cross section values that are consistent with observations of dwarf galaxy rotation curves, are in strong disagreement with each other. As long as this disagreement is unresolved, many of the additional constraints on the SIDM hypothesis described above, which also depend upon the validity of the N -body results or related analysis, will remain suspect. It is possible that either the N -body or the fluid model does not describe the system correctly. To remove such possibility, we directly compare the two methods in the simplest case, namely, the isolated spherically symmetric halo.
We will summarize the basis for the conducting fluid model in Section 2. In Section 3, we will describe our Monte Carlo numerical algorithm for the SIDM elastic scattering interaction in N -body simulations. We test those two methods against each other for isolated haloes in Section 4. The impact of our comparison results on the cosmological similarity solution by Ahn & Shapiro (2005) is discussed in Section 5. Finally, our results are summarized in Section 6.

Background: Gravothermal Collapse in Star Clusters
The conducting fluid model (also known as the gaseous model) was first developed to study gravothermal collapse in globular star cluster systems, and has been shown to be successful in those systems. Lynden-Bell & Eggleton (1980, hereafter LBE) proposed a thermal conduction formula for collisional, gravitationally-bound systems based on dimensional analysis, and derived an analytical solution that describes the self-similar collapse of a star cluster. That selfsimilar solution appears in the late stage of collapse in a Fokker-Planck calculation (Cohn 1980) and in N -body simulations (Makino 1996;Spurzem & Aarseth 1996;Baumgardt et al. 2003;Szell et al. 2005). When the time evolu-tion of a Plummer sphere is solved numerically by integrating the partial differential equations of the fluid model, the resulting collapse time agrees with that from other methods (Goodman 1987;Heggie & Ramamani 1989). SIDM haloes and globular clusters are both 'collisional,' self-gravitating systems, but the angular distribution and velocity dependence of the collisions are different in the two cases. Stars obey Rutherford scattering, which is dominated by smallangle scattering and small velocity encounters, σ ∝ v −4 , while the SIDM cross section we explore in this paper is isotropic and velocity independent. It is possible, however, that the SIDM interaction also obeys Rutherford scattering via a 'dark-photon' interaction (Ackerman et al. 2009;Feng et al. 2009). Gravitational Rutherford scattering between dark matter particles are negligible unless they are all 10 5 − 10 6 M mass black holes (Jin et al. 2005).

The Basic Equations
In this section, we review the conducting fluid model developed by LBE and BSI before we compare its results with our simulations. As in these papers, we shall here restrict the conducting model to the case of particle distributions that are spherically symmetric, isotropic in velocity dispersion, and quasi-static. The quasi-static approximation means that, while the fluid evolves thermally, it always satisfies hydrostatic equilibrium at each moment. For the problem at hand, it is a good approximation, because the collapse timescale is always much longer than the dynamical time. The conducting fluid model is not, in general, restricted to quasi-static systems, however (Bettwieser 1983;Ahn & Shapiro 2005). Deviation from isotropic velocity dispersion is also possible to consider, but it plays only a minor role for the problem of interest here. The Boltzmann equation, which is a partial differential equation in phase space, can be written as an infinite series of moment equations in position space by integrating over all velocities. When third moments are negligible (skewless), the series of moment equations can be closed by truncating at second order, as described in Ahn & Shapiro (2005). The equations which result in this case are what we shall here refer to as 'the fluid approximation.' The same equations would result from assuming a Gaussian velocity distribution (not necessarily isotropic or Maxwellian), an approximation called 'Gaussian closure' (Levermore 1996, and references therein). If velocity isotropy is imposed in addition, the fluid approximation gives the familiar Euler equation for an ideal gas with the ratio of specific heats, γ = 5/3. This fluid approximation can describe the structure formation of collisionless CDM, and accurately reproduce the CDM halo properties found in 3D N -body simulations when applied to the problem of spherical cosmological infall (Ahn & Shapiro 2005).
When collisions are important (e.g. either gravitational scattering between stars or DM self-interaction), thermal conduction, which is a third order moment, should not be neglected. However, accurate evaluation of third moments is only successful in the small Knudsen number regime, Kn 1. In this regime, the Navier-Stakes equation can be derived from the Boltzmann equation with the Fourier law of heat flux, where ρ(r, t) is the mass density, v(r, t) is the onedimensional velocity dispersion, σ is the scattering cross section per unit mass, λ ≡ 1/ρσ is the mean free path, and b ≡ 25 √ π/32 ≈ 1.38 is a constant 2 (Chapman-Enskog theory, e.g., Chapman & Cowling 1970;Lifshitz & Pitaevskii 1981). The local relaxation time is defined as tr(r, t) ≡ 1/(aρσv), where the constant a ≡ 16/π ≈ 2.26 describes the collision rate of particles that follow a Maxwellian distribution, defined in BSI.
In the other limit, Kn 1, Lynden-Bell & Eggleton (LBE) found that an empirical thermal conduction formula with λ replaced by the gravitational scale height (or Jeans length), H ≡ v 2 /4πGρ, explains the gravothermal catastrophe of star clusters very well, i.e., where C is an unknown constant of order unity. The scale height H characterizes the length scale that particles (or stars) orbit under the action of the gravitational force. BSI combined the two limiting forms of this thermal conduction into one as follows: The first term inside the brackets is the LBE formula, equation (2), which dominates in the large Knudsen number limit, λ H. In the other limit (λ H), the second term inside the brackets dominates, and the conduction converges to equation (1), instead. BSI's formula is an empirical interpolation between those two heat conductivity limits. BSI assumed C = b, but the exact value cannot be determined analytically. We determine the value of C by fitting the Monte Carlo N -body data in Section 4.1.1. The conducting fluid model is a set of moment equations closed empirically by this conductive heat flux L/4πr 2 , as follows: where M is the mass enclosed by radius r, and D/Dt is the Lagrangian derivative with respect to time. The first equation describes hydrostatic equilibrium. The second equation is the first law of thermodynamics, in the form relating energy transfer and entropy, 'dQ = T dS = dE + pdV .' When λ H, the fluid equations (4, 5) with the heat conduction equation (2) have a self-similar solution. In general, if a self-similar solution exists, if density is static as r → ∞ and if the evolution timescale ρc/ρc is proportional to the relaxation time at the centre, tr,c(t) ≡ tr(r = 0, t), then the central quantities evolve as, tr,c(t)/tr, for some constants α and t coll , by dimensional analysis (LBE). The exponents of 1 − t/t coll depend on the form of the relaxation time, therefore, different from those for star clusters. BSI obtained, by solving an eigenvalue problem of a system of ordinary differential equations.
For the non-self-similar time evolution considered in Sections 4.1.2, 4.1.3, and 4.2, we must integrate the time evolution of fluid variables ρ and v 2 numerically by alternative steps of the heat conduction and adiabatic relaxation to hydrostatic equilibrium, as described in BSI.

Scattering Algorithm
In this section, we will describe the Monte Carlo algorithm we implemented to model non-gravitational scattering of dark matter particles by other dark matter particles, within a pre-existing gravitational N -body method. Our scattering algorithm is similar to Kochanek & White (2000). Each particle can collide with one of its k nearest neighbors with a probability consistent with a given scattering cross section. For simplicity, we assume collisions are elastic, velocity independent, and isotropic in the centre of mass frame, but the Monte Carlo method can handle any differential cross section. We first outline the Monte Carlo N -body method and then explain, in detail, the algorithm that we have implemented in the parallel N -body code GADGET 1.1 (Springel, Yoshida, & White 2001), which uses the tree algorithm to calculate gravitational forces. Monte Carlo algorithms for particle-particle scattering (known as direct simulation Monte Carlo) have been used for more than thirty years to solve physics and engineering problems of collisional molecules, giving reasonable results (Bird 1994). For example, the results agree with an exact solution of the spatially homogeneous Boltzmann equation that describes the relaxation toward a Maxwellian distribution; they also agree with the Navier-Stokes equation solutions and experiments, including the thermal conductivity, in the small Kn regime (e.g., Nanbu 1984;Bird 1994;Gallis et al. 2004).
Consider N -body particles at positions xj and velocities vj with equal mass m. We discretize the distribution function f with, where W (x; r k ) is a spline kernel function of size r k , r kth j is the distance from particle j to its k ≈ 32nd nearest neighbor, and δ is the Dirac delta function. Our choice of kernel is often used in smoothed particle hydrodynamics, including GADGET. Our algorithm is identical to that of Kochanek & White if a top-hat kernel is used for W instead of a spline.
The result does not depend on the details of the kernel, however. We tested with k = 128 but did not see any difference.
The collision rate Γ for a particle at position x with velocity v to collide with this distribution f is, where σ is the scattering cross section per unit mass. Therefore the probability that an N -body particle 0 collides with particle j during a small timestep ∆t is, One can generate a random number and decide whether this collision happens and reorient velocities when they collide. This method is similar to a variant of direct simulation Monte Carlo called Nanbu's method (Nanbu 1980). His Monte Carlo algorithm, with the pairwise collision probability, equation (13), can be derived from the Boltzmann equation as described in his paper. Conversely, results of Nanbu's numerical method converge, mathematically, to the solution of the Boltzmann equation as the number of particles goes to infinity (Babovsky & Illner 1989). In Nanbu's method, only one particle is scattered per collision (only particle 0 but not j). The philosophy is that the N -body particles are samples chosen from real sets of microscopic particles, and those samples should collide with a smooth underlying distribution function, not necessarily with another sampled Nbody particle. However, then the energy and momentum are not conserved per collision. Moreover, the expectation value of the energy decreases systematically (Greengard & Reyna 1991). In our case, the error in the energy rises by 10 per cent quickly, so we decided to scatter N -body particles in pairs, not using Nanbu's method. Scattering in pairs is common in direct simulation Monte Carlo (e.g. Bird's method).
When particles are scattered in pairs, other particles j can scatter particle 0 during their timestep as well, but the scattering probability P0j, in equation (13), is similar to, but not exactly equal to Pj0, due to the difference in kernel sizes. Therefore, we symmetrize the scattering probability by taking the average scattering rate. Note that it is not trivial to generalize the pairwise scattering algorithm to simulations with unequal N -body particle masses, because P0j and Pj0 would then differ by a factor of their mass ratio; there is no reason to symmetrize two intrinsically different probabilities into single pairwise scattering probability.
In the following, we describe our algorithm in detail. Each particle, say particle 0, goes through the following steps, (i) to (iii), during its timestep ∆t0. Let particles 1, . . . , k be the k nearest neighbors of particle 0 (k = 32±2). The particle 0 collides with its neighbors with probabilities Pj0/2 (equation 13) during its timestep. The factor of two is the symmetrization factor that corrects the double counting of pairs. A particle j would also scatter particle 0 during its timestep, which results in a symmetrized scattering rate. Imagine a probability space [0, 1], with disjoint subsets Ij ≡ [ j−1 l=1 P l0 /2, j l=1 P l0 /2) that represent scattering events between particles 0 and j. We neglect the possibility of multiple scattering in the given timestep. Particle 0 collides with at most one of its neighbors. We restrict the timestep so that it is small enough that this approximation is good enough (see equation 15 below). We generate a uni-form random number x in [0, 1], and scatter particles 0 and j if x falls in a segment Ij, as described below.
(i) In the large Kn regime, most of the particles do not collide with another particle in a given timestep. Therefore, we can reduce the computation by estimating the rough scattering probability first, and compute the accurate probability P0j only if necessary. First, we calculated an upper bound to the scattering probability, whereρ is the approximate density calculated from r kth 0 viã ρ = km/ 4 3 π(r kth 0 ) 3 , and vmax is the maximum speed of all the particles. If the generated random number x is larger than P , this means that x is not in any segment Ij, therefore, the particle 0 does not collide during this timestep.
(ii) If the possibility of collision was not rejected in step (i), we calculate the pairwise scattering probability Pj0 and determine which neighbor the particle 0 collides with. The index j of the collision partner is the smallest integer that satisfies x j l=1 P0j, i.e. x ∈ Ij, if such j exists (otherwise, the particle does not collide with any neighbors).
(iii) For particle pairs that collide, we reorient their velocities randomly, assuming an elastic scattering which is isotropic in the centre of mass frame. Isotropic random directions can be generated by one square root operation, without using trigonometric functions (e.g., Vesely 2001). The velocities are updated in the kick phase of the leap-frog time integration in the GADGET 1.1 gravitational N -body method. At that time, we also update the centre-of-mass velocities of the nodes around the scattered particles in the oct-tree, used for the gravity calculation. This is because the centreof-mass velocities can be changed drastically by scattering.
We allowed at most one scattering per timestep per particle. In order to suppress the error due to possible multiple scattering, we restrict the timestep so that, This restriction makes the Monte Carlo method computationally costly in the small Kn regime, because then the timestep becomes much smaller than the dynamical timeof the order of the timestep in collisionless N -body simulations. This timestep problem may seem to be avoided by performing multiple scatterings per dynamical timestep, but there is another limit when the Knudsen number is small. The distances to kth neighbors must be smaller than the mean free path λ = 1/(ρσ), Otherwise, particles more than a mean free path away would be allowed to collide and, thereby, make the heat transfer larger than it should be in the diffusion limit. If one tries to avoid this by choosing a kernel size smaller than the mean particle separation length, then the particles simply freely stream beyond the mean free path, which is again incorrect. The only way to overcome this problem is to increase the number of the particles in inverse proportion to the volume within the mean free path, N particles ∝ λ −3 = (ρσ) 3 , which increases very rapidly during core collapse as ρ increases.
Conditions (15) and (16) prevent us from running simulations with small mean free path.

Dimensions and Units
We use the following characteristic scales in the rest of this paper. For the initial BSI self-similar profile (Sections 4.1.1, 4.2) and the Plummer model (Section 4.1.2), we describe densities and velocities in units of central quantities, ρc(t) ≡ ρ(0, t) and v 2 c (t) ≡ v 2 (0, t); ρ(r, t) and v(r, t) denote the density and the one dimensional velocity dispersion, respectively, in spherical symmetry. We use core radius rc(t) ≡ vc/ √ 4πGρc as a standard length scale. For Hernquist and NFW profiles (Section 4.1.3), which have singularities at the centre, we use the scale radius rs and density ρ0 that appear in the density profiles as units, instead. We use v0 ≡ rs √ 4πGρ0 as the velocity scale, which is similar to the definition of the core radius above.
The relaxation times at the centre tr,c(t) ≡ tr(0, t) for BSI and Plummer model, or tr,0 ≡ 1/(aρ0σv0) for NFW and Hernquist profiles, are used as unit time-scales.
We express cross sections in a dimensionless way aŝ σc ≡ ρcσrc andσ0 ≡ ρ0σrs, where σ is the scattering cross section per unit mass. The inverses are the Knudsen numbers, the mean free path in the unit of system size (rc or rs). The evolution of the system is characterized by the Knudsen number Kn only, not depending on the overall physical scale.

Simulation Setup
We generate the initial conditions for N -body particle positions and velocities randomly from the distribution functions using the rejection method (Aarseth et al. 1974). 3 The distribution function of the BSI's self-similar solution and the NFW profile are calculated numerically by Eddington's formula (Binney & Tremaine 1987). The distribution functions of the Plummer model and the Hernquist model (Hernquist 1990) have known analytical forms. We set the initial centreof-mass velocity of the system of N -body particles to zero by an overall boost.
We truncate the initial profile at some radius r f , and put a simple reflecting boundary, which flips the direction of the radial velocity if particles are moving outward passed the reflecting boundary. We use r f = 600 rc for the BSI self-similar profile, r f = 58.5 rc for the Plummer model, and r f = 100 rs for the Hernquist and NFW profiles. The density at r f is smaller than 2 × 10 −7 ρc, and the heat flux (L/4πr 2 in Section 2) at r f , calculated from the equation for the conducting fluid model, is smaller than 0.02 per cent of its maximum value. Nevertheless, the collapse time is sometimes sensitive to the position of this reflecting boundary. The collapse was slower by a factor of two when we first used truncation radius r f = 58.5 rc for the BSI profile, even though this r f is large enough to satisfy Antonov's criterion for gravothermal catastrophe (Endoh et al. 1997). We also tested r f = 300 rc, for the BSI profile, and the collapse time changed by only 3 per cent compared to r f = 600 rc. Hence, our choice of r f = 600 rc is large enough to bring about a converged solution to high accuracy.
We use two timestep criteria in addition to equation is the maximum number of particles inside the core, r f is the radius of the reflecting boundary, is the initial Plummer equivalent gravitational softening length, which is reset sometimes overtime as central density ρc(t) increases; i.e. ∝ rc(t) ∝ v 2 c / √ ρc. The constant C is the LBE prefactor used for the conducting fluid model.
(15). The timestep must also satisfy ∆t ηvvc(0)/a, and ∆t ηG/ √ Gρ, where vc(0) is the initial one-dimensional velocity dispersion, a is the local acceleration, andρ is the local density calculated from k = 32 nearest neighbors (see below equation 14). We choose the dimensionless parameters to be ηv = 0.02 and ηG = 0.005. With this choice, the initial conditions are static for several dynamical times when scattering is turned off. Energy conservation is satisfied to better than 1 per cent in all runs.
The number of particles, gravitational softening length and other parameters are summarized in Table 1. For cases BSI and P, we reset the gravitational softening length to 0.1 rc(t) every time the central density increases by a factor of 10, to avoid the numerical effect of softening on the density profile.
We tested that the numerical scattering rate is correct, by counting the number of scatterings in the simulation of a non-singular isothermal sphere 4 and comparing it to the analytical rate. The scattering rates agree within 3 per cent when 64 3 particles are used and the isothermal sphere is truncated at 58.5 rc. The difference is due to the fluctuation in the randomly generated initial condition, not to the Poisson fluctuation in the number of scatterings.

Analysis Methods
We calculate the central quantities for each snapshot from N -body particles within a sphere of radius rc around the density-weighted centre of mass (Casertano & Hut 1985), assuming an isothermal sphere inside rc. Namely, we find a consistent rc iteratively that satisfies, where ρc is the central density estimated from M (rc), the mass inside rc, and vc is the velocity dispersion inside rc. The value 1.10 in equation (18) is the ratio of the central density to the average density inside rc for the non-singular isothermal sphere. In this way, we can use as many particles as possible without systematic error for the calculation of the central quantities. The velocity dispersion inside the core quickly becomes isothermal due to collisions.
We calculated the smoothed density and velocity dispersion field at each point in space using an adaptive kernel, similar to that used in the well-known SPH method, with a Gaussian kernel whose size is r32 -the distance to the 32nd nearest neighbor. These smoothed-particle density and velocity dispersion fields were then averaged over spherical surfaces on different radii r to give the radial profiles of these quantities. In practice, this amounts to summing the spherically-averaged Gaussian kernels evaluated at each radius r (e.g., Reed et al. 2005).

RESULTS
We compare our Monte Carlo N -body simulations with the conducting fluid model, first in the large Kn regime for various initial conditions, and then in the transitional regime Kn ∼ 1.

Self-Similar Gravothermal Collapse Solution
In this section, our Monte Carlo N -body simulation for large Kn,σc(0) = Kn −1 = 0.067, is compared with the BSI selfsimilar solution. Fig. 1 shows the evolution of the density and velocity dispersion profiles. The right panels, plotted in self-similar variables, show that, when the N -body particles are initialized according to a given time-slice of the self-similar solution, they indeed evolve self-similarly in the Monte Carlo N -body simulation, thereafter. Fig. 2 shows that our Monte Carlo N -body simulation is in excellent agreement with the self-similar solution (equations 6-10), by adjusting the value of one parameter, C, in the heat conduction equation (2). We determine the collapse time t coll = 385 tr,c(0) by fitting the time evolution of the relaxation time data (Fig. 2, left-bottom panel) by the linear function, equation (8). The prefactor of the thermal flux C = 0.75 follows from equation (10). The best-fitting power law index of v 2 c ∝ ρ (α−2)/α c (equations 6 and 7) gives a value α = 2.22 (Fig. 2, right-bottom), which agrees reasonably well with the value of the self-similar solution, 2.19 (equation 9).
To test convergence, we simulated the self-similar solution with three different numbers of particles, 2×32 3 , 4×32 3 , and 16×32 3 , which give the initial number of particles inside the core Nc(0) = 227, 376 and 1527, respectively. At the time of our resolution test, we used reflecting radius r f = 58.5 rc. Numerical errors due to finite number of particles should only depend on the N -body particle density, independent of    the choice of boundary condition. As a result, it was not necessary to run additional convergence test for the case with r f = 600 rc, used in our final runs. Other parameters are the same as run BSI. In Fig. 3, we plot the evolution of density and number of particles inside the core as a function of time.
The three smooth curves in the left panel are the same analytic self-similar solution (equation 6) with t coll = 705 tr,c(0). The results for the runs with N = 4 × 32 3 and N = 16 × 32 3 agree very well, while those for N = 2 × 32 3 deviate from the other two runs systematically when Nc 100. Our run BSI contains 470 particles inside the core at t = 0 (Table 1), which is better than the converged N = 4 × 32 3 run here.
Other runs in the following sections have better resolution inside the core, because of the smaller fraction of particles for r > rc, resulting from the steeper decline in their density profiles.

Plummer Model
We compare the Monte Carlo N -body simulation with the conducting fluid model when the initial condition is the Plummer model in the large Kn regime. The Plummer model, a standard initial condition used to study gravothermal instabilities, has a spherical mass distribution given by, In this case, the characteristic scales defined in Section 3.2 can be shown to be, vc = (2GMT /a pl ) 1/2 /12 and rc = a pl /(3 √ 2). We evolve this system according to the conducting fluid model in the large Kn limit with the quasi-static approximation described in Section 2, equations (2), (4), and (5). The N -body simulation is also performed in the large Kn regime,σc(0) = 0.013. The time evolution for this quasistatic system is independent of the actual value of σ for the large Kn regime, if the time for solutions with different σ is expressed in units of the relaxation time.
We plot the time evolution in Fig. 4. The fluid model agrees with the N -body results reasonably well if the coefficient C is assigned a value of 0.80. This is in reasonably good agreement with the value of C = 0.75, found for the self-similar solution. The logarithmic slope (right-top panel) has a plateau at about −α, which is the asymptotic slope of the self-similar profile. This implies that the inner part is converging to the self-similar solution, with an asymptotic logarithmic slope −α, which was well known in the gravothermal collapse of star clusters.

NFW & Hernquist Profiles
We also considered the case of an initial condition with spherical mass distribution given by the NFW profile (Navarro et al. 1997) to test the consistency of the Monte Carlo N -body simulation and the conducting fluid model with each other for the typical halo profile seen in cosmological N -body simulations. The gravothermal collapse timescale is interesting as a way to put constraints on the SIDM cross section, because the result of gravothermal collapse is a divergent profile with ρ ∝ r −2.2 , which is not seen in observations. We note that gravothermal collapse is prevented by cosmological infall or major mergers, but when a halo decouples from cosmological growth, the halo could experience gravothermal collapse. The NFW profile is given by, In addition, we perform a test with initial conditions based upon a Hernquist profile (Hernquist 1990), which has the same inner profile, ρ ∝ r −1 , and is sometimes used to approximate the NFW profile. As we will demonstrate below, the gravothermal collapses for the NFW and Hernquist profiles are significantly different, and it is therefore, dangerous to use the results for the Hernquist profile (Kochanek & White 2000) to describe realistic haloes. As for the Plummer model, our conducting fluid model calculations adopt the large Kn limit, and the N -body simulations use a small, but nonzero, σ (cf. Table 1). Figs. 5 and 6 show that the fluid model agrees reasonably with N -body results. In particular, the density profiles in the two cases evolve through a sequence which is close to the N -body results, although the agreement is better if the following adjustments are made to the time axis of the fluid model solutions (6). For NFW, N -body and fluid model central density agree best when compared at times which differ by a small shift equal to 100 tr,0, or 20 per cent of the collapse time. For the Hernquist initial profile, the agreement is best if the time axis is scaled by a factor of order of unity which is equivalent to adjusting the value of C in the conductivity from C = 0.75 (found in the BSI run) to 0.9, which also differs by only 20 per cent. For both NFW and Hernquist initial conditions, density profiles agree very well when the central densities are equal, but the values of the central densities differ by about 20 per cent at the maximum core expansion.
We note that the initial NFW profile evolves very differently from the initial Hernquist profile, for the same parameters ρ0 and rs. The NFW run has a central density about three times smaller at maximum core expansion, and a collapse time about four times longer than for the H run. This is because the NFW profile has larger heat flux at r rs due to its larger density there, which heats and expands the central mass more than does the Hernquist profile. . We use C = 0.75 for NFW and C = 0.9 for Hernquist profile. In the left panel, the origin for the fluid curve is shifted to the right along the time axis by 100 t r,0 to show the agreement in the gravothermal collapse phase. Minimum density of the fluid is ±20 per cent different from N -body. Error bars show ∆ρc/ρc = 2/ √ Nc for every 10 data points.
Our simulation results are qualitatively similar to those of Kochanek & White (2000), but the time evolutions differ by a factor of two. We do not know the reason for this difference. We show how our units convert to those in Kochanek & White (2000) as follows: ρ KW 0 = 2πρ0, t KW r,c = 1.7 × 2 √ 2atr,0 = 11 tr,0 andσ KW DM = 2πσ0. Their simulation forσ KW DM = 1, which has the same cross section as ours, reached minimum density at t ≈ t KW r,c = 11tr,0 and collapsed gravothermally to ρc = 2ρ KW 0 = 13 ρ0 at t ≈ 3.2t KW r,c = 35 tr,0, while our simulation reached those densities at 20 tr,0 and 70 tr,0, respectively (Fig. 6). Their evolution is about twice as fast as in our simulation. The origin of this discrepancy is unknown.

The Source of Difference
Our Monte Carlo N -body simulations agree with the solutions of the conducting fluid model in the self-similar gravothermal collapse phase, but in general have about a 20 per cent difference in the central density and the collapse rate. This modest disagreement is not surprising in view of the fact that the conducting fluid model is not an exact theory derived from first principals. Heggie & Stevenson (1988) compared the thermal conductivity of the conducting fluid model with that calculated from the orbit-averaged Fokker-Planck equation for star clusters with several profiles, including polytropes and lowered Maxwellians, evaluated at the centres. They find that the coefficients of conductivity C varied from one profile to another by factors of 2 or 3. The overall collapse rates are not that different, probably because the profiles quickly converge to the self-similar solution around the centre. Indeed, for star clusters, the value of C = 0.88, which makes the fluid model match the asymptotic collapse rate in the isotropic Fokker-Planck calculation (Cohn 1980), also gives the correct collapse time of the Plummer model, 15.4 half-mass relaxation times (Goodman 1987;Heggie & Ramamani 1989). This means that C needs not to be different for the cases of self-similar collapse solution and the collapse of Plummer profile.
One possible distinction between the 1D conducting fluid model and the 3D N -body/Monte Carlo simulations is that the former assumes that particle velocity distribu-q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 2 v 2 v c 2 (0) BSI q radial tangential q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.1 1 10 Plummer q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 0.1 v 2 v 0 2 NFW q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q tions are isotropic in the frame of the bulk flow while the latter do not. Anisotropy in velocity dispersion affects the collapse time (Bettwieser 1983;Louis 1990). Fokker-Planck calculations with anisotropy show that the collapse time of the star clusters initialized by the Plummer model is 20 per cent larger than for the isotropic case (Takahashi 1995), and agrees with N -body simulation (Khalisi et al. 2007). In Fig. 7, we plot the radial (v 2 r ) and tangential [v 2 ⊥ = (v 2 θ + v 2 φ )/2] velocity dispersions of our simulation when the central density has increased by a factor of 10. We do not see anisotropies near the centre, but we do at r 5rc for runs P, NFW and H. When particles scatter from the centre to large radii, their nearly radial orbits bring anisotropy to the initially isotropic velocity distribution at large radii. If the original unscattered particles at large radii are less numerous, because the density profile is steeper, the anisotropy that results from scattering particles from small to large radii will be relatively larger. This is why anisotropy is larger if the logarithmic slope of the density profile is steeper. Anisotropy may therefore play some role in the collapse rate of run P and H.
In short, the conducting fluid model has been shown to describe the gravothermal collapse relatively accurately, but we cannot expect very high precision in general. Therefore, the 20 per cent match in minimum core density of runs NFW and H, and 20 per cent match in the gravothermal collapse time for runs P, NFW and H run is a very reasonable agreement.

Transitional Regime
When Kn becomes comparable to, or smaller than, one, the self-similar collapse solution of BSI no longer applies, because of the presence of the second term in the heat conduction equation (3). The time evolution in units of the relaxation time becomes slower the smaller Kn gets, because smaller mean free path suppresses particle transport more, which results in smaller heat conduction (BSI). To test the heat conduction equation (3) (0) is the ratio of the mean free path to the kernel size at the centre at t = 0, and t 10 is the time for the density to increase by a factor of 10. Other parameters are the same as for the BSI run.
body simulations with those of the conducting fluid model, derived numerically when no self-similar solution exists, with initial Knudsen numbers Kn −1 (0) =σc(0) = 0.25, 0.5, 0.75 and 1.0. The initial condition is the same as in Section 4.1.1: the particle distribution of the BSI self-similar collapse solution. See Table 2 for the summary. N -body simulation with largerσc becomes very difficult because of the mean free path requirement equation (16). The initial ratio of mean free path to the kernel size, λ/r32 ∝ ρ −2/3 , at the centre is listed in the table. Ourσc(0) = 1.0 run violates λ > r32 at ρc(t) ∼ 8. The number of particles required to follow the collapse scales asσ 3 c beyond this density or cross section. For the fluid model, the prefactor C is chosen to be 0.75 to make an agreement at smallσc (Section 4.1.1). We ran the fluid code with initial dimensionless cross sectionσ(0)/ √ b = 0, 0.25, 0.5, 0.75, 1.0, 1.5 and 2.0. The time evolution of the fluid model in units of initial relaxation time depends only on the combinationσ/ √ b, which can be seen from the heat conduction equation (3).
We plot the evolution of the central density for both N -body and fluid model results in Fig. 8. In units of the relaxation time, the collapse is slower for larger cross section. Fluid model results forσc(0)/ √ b = 0.5, 1.0, 1.5 and 2.0 evolve similarly to N -body simulation results forσc(0) = 0.25, 0.5, 0.75 and 1.0, respectively. This suggest that the effective value of the coefficient is b = 0.25 in the transitional regime, Kn 1.
To show the agreement between our N -body results and the fluid model solution with b = 0.25, we plotted the normalized collapse rate as a function of cross section in Fig. 9, where t10 is the time for the density to increase by a factor of 10, and the fiducial time scale t * 10 is the time at which the self-similar solution withσc(0) = 1 has a density increase by a factor of 10. The collapse rate is proportional to the cross section in the large Kn regime (smallσ), but deviates from the linear relation in the transitional regimê σ 0.5, as predicted by BSI. The collapse rate should reach some maximum at some cross section and then decrease as t −1 10 ∝σ −1 c asσc → ∞. Furthermore, since the value of b can be calculated from first principals in the small Kn regime, b should converge to the Chapman-Enskog value (b=1.38 in Section 2) as Kn → 0. However, due to the numerical limit in equation (16), we cannot go into the small Kn regime to see the convergence to the Chapman-Enskog theory or the turn over of the collapse rate. Our N -body results are consistent with a constant b = 0.25 in the range we are able to simulate.  Until now, we have limited our attention to the gravothermal relaxation of isolated SIDM haloes. However, the haloes that form in a cosmological context are not isolated but rather build up over time from the nonlinear growth of small-amplitude initial density perturbations in the expanding cosmological background universe, which results in a continuous infall of additional mass. To model this cosmological formation and evolution of SIDM haloes, A&S added the BSI heat conduction term to the fully timedependent conservation equations of the fluid approximation in 1D, spherical symmetry, as described in A&S. They used these equations to derive an analytical solution for the so-called 'secondary infall' problem, in an Einstein-de Sitter universe, with initial perturbation given by the spherical overdensity profile δM/M ∝M −1/6 , whereM (r) is the unperturbed mass in a sphere of radius r at the mean density of the universe, and δM (r) = M (r) −M (r) is the mass perturbation inside the radius. The solution with this initial overdensity is self-similar, that is, the solution is timeindependent if radius and density are measured in units of the time-varying turn-around radius rta and background critical density ρ b , respectively. The family of similarity solutions is parametrized by the dimensionless cross section, where rvir is the halo virial radius, the radius at which an accretion shock occurs. The density profile has a core whose density and size depend upon the value of Q. The heat conduction equation (3) has a dependence on σ = H/λ given by, This function takes its maximum value √ abC/2 atσ = √ a −1 b C −1 . Compared to previously assumed values, C = b = 1.0, our calibrated heat conduction has a maximum that is smaller by a factor of √ bC ≈ 0.4, which occurs at a value ofσ that is smaller by a factor of √ bC −1 = 0.58. There is a minimum central density and largest core size which occur approximately when the dimensionless cross section at the centre,σc, defined in Section 3.2, maximizes the heat flux with the value,σc = √ a −1 b C −1 . The solution with thatσc is defined as the maximally-relaxed solution, and the corresponding cross section Q is denoted by Q th (A&S). This maximally-relaxed halo has a density profile almost identical to the empirical Burkert profile (Burkert 1995), which fits the observed rotation curves of dwarf and LSB galaxies well.
We used the same numerical code as A&S to solve for the cosmological similarity solutions with our calibrated prefactors, and to compare with their original solutions. In Fig. 10, we plot the maximally-relaxed density profile (left panel), and the dependence of the central density on dimensionless cross sectionσc (right panel). For prefactors C = b = 1.0 adopted by A&S, the solution is maximally relaxed for Q th = 7.35 × 10 −4 , with central density ρc = 1.17×10 4 ρ b . For C = 0.75 and b = 0.25, the maximallyrelaxed solution shifts to Q th = 3.95 × 10 −4 , with a central density ρc = 1.50 × 10 4 ρ b . 5 This central density is about 30 per cent larger than that of the original solution. For σc 1, a cross section that is 1.33 times as large is required to achieve the same central density, due to the change of coefficient C from 1.0 to 0.75. 5 The value ofσc which defined the maximally-relaxed solution, σc = √ a −1 bC −1 , changes from 0.666 to 0.384 if the conductivity parameters change their values in A&S to our new values here. However, the actual minimum densities occur at slightly higher σ. For C = b = 1.0, central density takes its minimum value ρc = 1.15 × 10 4 ρ b atσc = 0.85, or Q = 9.78 × 10 −4 . For C = 0.75 and b = 0.25, the minimum density ρc = 1.46 × 10 4 ρ b is found at σc = 0.55, or Q = 3.95 × 10 −4 . A&S estimated the cross section values that brought dwarf galaxy rotation curves into agreement with the cosmological self-similar halo profiles. In order to translate a given value of Q th into a cross section value σ, the typical formation time for haloes that host dwarf galaxies in the ΛCDM model of structure formation was used to relate rvir and ρ b in the definition of Q, as described in A&S. Using the same argument, our corrected Q th value corresponds to a cross section 117 cm 2 g −1 , which is smaller than the original value 218 cm 2 g −1 but still much larger than the values, 0.5 − 5 cm 2 g −1 , found in Monte Carlo N -body simulations that produce cored profiles for galactic haloes formed from cosmological initial conditions (Davé et al. 2001).
We will discuss this apparent discrepancy between the values of the SIDM cross section which are best able to match the observed rotation curves of dwarf and LSB galaxies, from the A&S similarity solutions and cosmological Nbody/Monte Carlo simulations, respectively, in a separate paper. Here, we have demonstrated that this discrepancy is not likely to result from the break down of the conducting fluid model. We have found good agreement, in fact, between the Monte Carlo N -body simulations and the conducting fluid model for the thermal relaxation and gravothermal collapse of isolated haloes (of fixed mass), at least with regard to the long mean-free-path regime and the transitional regime in which the mean free path is comparable to the system size. These are the regimes of greatest relevance to the SIDM halo problem. We have also shown that, when we use the agreement between our Monte Carlo N -body results and the conducting fluid model solutions for isolated haloes to calculate the dimensionless parameters on which the heat conduction depends, the impact of the modified parameters on the A&S similarity solution is relatively small. We must, therefore, seek a different explanation for the different values of σ required by the cosmological N -body/Monte Carlo and the similarity solutions of A&S, to produce the same observed degree of relaxation of the SIDM halo density profiles.
As we shall discuss in our next paper in this series, the essential difference between the self-similar solution of A&S and the cosmological N -body simulations is that the self-similar halo is continuously heated by the accretion shock, while the inner part of a more realistic halo, on the galaxy scale, is eventually unaffected by infall when infall rate slows after the halo forms. We will compare cosmological Monte Carlo N -body simulations with the conducting fluid model with non-self-similar infall in our next paper.

CONCLUSION AND SUMMARY
The ability of the SIDM hypothesis to resolve the cusp-core problem of collisionless CDM haloes on the galaxy scale and the cross section values required to accomplish this remain uncertain, as long as the previous conclusions drawn from N -body simulation and the conducting fluid model differ so strongly. To rectify this situation, we have developed a new Monte Carlo N -body code of our own, based on the pre-existing GADGET 1.1 N -body code, and applied it to compare the N -body simulations with the conducting fluid model for isolated, spherically symmetric self-gravitating SIDM haloes. The collisions were assumed to be velocity independent, elastic, and isotropic. Our results include the following: (i) Our Monte Carlo N -body simulations are in very good agreement with the analytical self-similar gravothermal collapse solution of BSI, when the coefficient of thermal conduction is set to C = 0.75 (Section 4.1.1); the density and velocity dispersion profiles evolve self-similarly; the central density and velocity dispersion follow the self-similar time evolution formulae with the predicted constant α = 2.19; and, the time to collapse is always proportional to the central relaxation time at that time (equation 6-10).
(ii) The conducting fluid model agrees with Monte Carlo N -body simulations reasonably well for different initial conditions: Plummer model, Hernquist profile, and NFW profile, in the large-Kn regime. The collapse time and the central density at maximum core expansion agree within 20 per cent. The shape of the density profile and the central density evolution as a function of time during gravothermal collapse agree very well.
(iii) We also showed that the collapse time becomes longer in units of relaxation time as the system transitions from the large-to the small-Kn regime, as predicted by the conducting fluid model. The N -body results agree with the conducting fluid model for Kn 1, orσ 1, with the prefactor b = 0.25. However, this prefactor is more than five times smaller than the Chapman-Enskog value, valid asymptotically in the small Kn limit. The conducting fluid model must be further calibrated against N -body simulations if it is used beyond the transitional regime, forσ 1.
(iv) Our calibration of the prefactors C and b does not change the cosmological similarity solutions of Ahn & Shapiro (2005) significantly. The cross section that gives the minimum central density on the dwarf galaxy scale is altered from 220 cm 2 g −1 to 117 cm 2 g −1 , but this is still much larger than the values that cosmological Monte Carlo Nbody simulations used to make cored SIDM haloes (σ ∼ 0.5 − 5 cm 2 g −1 ). We will investigate this problem in our subsequent paper. Our results here suggest that this apparent discrepancy is not the result of a breakdown of either the conducting fluid model or the Monte Carlo scattering algorithm in the N -body simulations. As we shall show in a com-panion paper, in fact, the discrepancy results, instead, from the gradual departure of halo evolution from self-similarity as the infall rate during cosmological structure formation drops below the self-similar rate at late times after halo formation.